stanEx2
ex5 regression
ex5.stan
data{
int N;
int K;
vector[N] y;
matrix[N,K] X;
}
parameters{
vector[K] b;
real<lower=0> s;
}
model{
vector[N] m=X*b;
y~normal(m,s);
}
generated quantities{
vector[N] y1;
vector[N] m1=X*b;
for(i in 1:N){
y1[i]=normal_rng(m1[i],s);
}
}normal linear models
n=30
mdl=cmdstan_model('./ex5.stan')single regression
tb=tibble(x=runif(n,0,9),
y=rnorm(n,x,1))
f=formula(y~x)
par(mfrow=c(1,1))
plot(tb$x,tb$y)
qplot(data=tb,x,y)
X=model.matrix(f,tb)
data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X)
mle=mdl$optimize(data=data)## Initial log joint probability = -646.428
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 15 -13.4147 0.00105322 0.000473195 1 1 23
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -13.41
## b[1] -0.11
## b[2] 1.00
## s 0.95
## y1[1] 5.19
## y1[2] 9.59
## y1[3] 4.63
## y1[4] 5.37
## y1[5] 4.62
## y1[6] 6.14
## y1[7] 8.18
## y1[8] 6.52
## y1[9] 3.68
## y1[10] 1.32
## y1[11] -0.74
## y1[12] 5.57
## y1[13] 7.96
## y1[14] 2.44
## y1[15] 1.89
## y1[16] 0.68
## y1[17] 7.27
## y1[18] 4.85
## y1[19] 5.32
## y1[20] 4.03
## y1[21] 5.97
## y1[22] 1.60
## y1[23] 5.44
## y1[24] 6.08
## y1[25] 0.90
## y1[26] 6.35
## y1[27] 3.06
## y1[28] 2.60
## y1[29] 7.14
## y1[30] 3.89
## m1[1] 5.92
## m1[2] 8.31
## m1[3] 5.58
## m1[4] 5.05
## m1[5] 4.43
## m1[6] 6.15
## m1[7] 8.33
## m1[8] 7.61
## m1[9] 4.53
## m1[10] 1.00
## m1[11] 0.07
## m1[12] 6.96
## m1[13] 7.58
## m1[14] 1.84
## m1[15] 1.62
## m1[16] 0.48
## m1[17] 7.64
## m1[18] 4.03
## m1[19] 4.00
## m1[20] 4.95
## m1[21] 5.15
## m1[22] 1.18
## m1[23] 8.66
## m1[24] 7.98
## m1[25] 0.89
## m1[26] 7.73
## m1[27] 3.24
## m1[28] 3.40
## m1[29] 8.27
## m1[30] 2.88mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.4 seconds.seeMCMC(mcmc,exc='m',ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -15.00 -14.65 1.29 1.00 -17.71 -13.64 1.01 644 531
## b[1] -0.12 -0.11 0.40 0.39 -0.81 0.50 1.01 382 466
## b[2] 1.00 1.00 0.07 0.07 0.89 1.11 1.01 453 641
## s 1.02 1.00 0.14 0.13 0.81 1.30 1.01 1115 947
## y1[1] 5.88 5.87 1.02 1.02 4.18 7.58 1.00 2019 1696
## y1[2] 8.33 8.31 1.06 1.05 6.59 10.10 1.00 1796 2104
## y1[3] 5.59 5.57 1.05 1.05 3.88 7.37 1.00 2037 1927
## y1[4] 5.05 5.06 1.07 0.98 3.28 6.80 1.00 1978 1869
## y1[5] 4.43 4.44 1.07 1.01 2.67 6.14 1.00 2033 1970
## y1[6] 6.13 6.12 1.07 1.04 4.38 7.88 1.00 2100 1963
## y1[7] 8.31 8.32 1.06 1.05 6.57 10.09 1.00 1975 1933
## y1[8] 7.61 7.62 1.08 1.06 5.85 9.37 1.00 1946 1929
## y1[9] 4.51 4.52 1.05 1.01 2.83 6.24 1.00 2119 1952
## y1[10] 1.00 0.99 1.06 1.04 -0.76 2.72 1.00 1534 1968
## y1[11] 0.08 0.07 1.14 1.10 -1.86 1.92 1.00 1367 1825
## y1[12] 6.96 6.96 1.03 1.04 5.29 8.68 1.00 1788 1856
## y1[13] 7.61 7.61 1.06 1.05 5.84 9.31 1.00 2169 1802
## y1[14] 1.85 1.85 1.09 1.06 0.05 3.69 1.00 1614 1777
## y1[15] 1.60 1.61 1.11 1.04 -0.21 3.46 1.00 1744 1704
## y1[16] 0.49 0.48 1.11 1.01 -1.36 2.28 1.00 1433 1965
## y1[17] 7.65 7.64 1.08 1.06 5.92 9.47 1.00 1838 1772
## y1[18] 3.99 4.00 1.03 0.98 2.30 5.68 1.00 1860 1774
## y1[19] 4.05 4.05 1.06 1.02 2.29 5.79 1.00 1895 1933
## y1[20] 4.96 4.94 1.07 1.06 3.25 6.73 1.00 2030 1950
## y1[21] 5.18 5.19 1.05 0.99 3.48 6.90 1.00 1945 1849
## y1[22] 1.17 1.19 1.10 1.04 -0.66 2.91 1.00 1543 1591
## y1[23] 8.64 8.62 1.05 1.07 6.92 10.33 1.00 1860 1646
## y1[24] 8.00 7.99 1.06 1.05 6.30 9.80 1.00 1757 1760
## y1[25] 0.91 0.92 1.08 1.04 -0.84 2.66 1.00 1641 1352
## y1[26] 7.71 7.73 1.06 1.01 5.94 9.40 1.00 2057 1926
## y1[27] 3.25 3.28 1.04 1.02 1.50 4.91 1.00 2134 1908
## y1[28] 3.38 3.44 1.04 1.00 1.64 4.99 1.00 1961 1880
## y1[29] 8.26 8.27 1.05 1.07 6.51 9.98 1.00 2169 1930
## y1[30] 2.85 2.87 1.04 1.03 1.15 4.50 1.00 1741 1812
## m1[1] 5.92 5.91 0.20 0.20 5.61 6.27 1.00 2101 1646
## m1[2] 8.31 8.31 0.30 0.29 7.83 8.81 1.00 1266 1156
## m1[3] 5.58 5.57 0.19 0.19 5.27 5.90 1.00 1808 1460
## m1[4] 5.05 5.05 0.19 0.19 4.75 5.37 1.00 1283 1231
## m1[5] 4.42 4.42 0.19 0.19 4.12 4.74 1.00 837 1012
## m1[6] 6.15 6.14 0.20 0.20 5.83 6.50 1.00 2195 1601
## m1[7] 8.33 8.33 0.30 0.29 7.85 8.83 1.00 1258 1149
## m1[8] 7.61 7.60 0.26 0.25 7.19 8.05 1.00 1657 994
## m1[9] 4.53 4.53 0.19 0.18 4.23 4.85 1.00 903 1033
## m1[10] 0.99 1.01 0.33 0.32 0.42 1.50 1.01 387 478
## m1[11] 0.07 0.08 0.39 0.38 -0.61 0.67 1.01 383 466
## m1[12] 6.96 6.96 0.23 0.23 6.59 7.35 1.00 2045 1366
## m1[13] 7.58 7.58 0.26 0.25 7.17 8.02 1.00 1680 994
## m1[14] 1.84 1.85 0.28 0.28 1.34 2.27 1.01 400 497
## m1[15] 1.61 1.63 0.30 0.29 1.10 2.06 1.01 398 496
## m1[16] 0.47 0.49 0.36 0.36 -0.16 1.03 1.01 384 472
## m1[17] 7.64 7.64 0.26 0.25 7.22 8.08 1.00 1634 994
## m1[18] 4.03 4.02 0.20 0.19 3.72 4.35 1.01 669 921
## m1[19] 4.00 4.00 0.20 0.19 3.69 4.32 1.01 661 895
## m1[20] 4.95 4.95 0.19 0.19 4.65 5.27 1.00 1200 1200
## m1[21] 5.15 5.15 0.19 0.19 4.85 5.47 1.00 1371 1232
## m1[22] 1.18 1.19 0.32 0.32 0.63 1.67 1.01 390 478
## m1[23] 8.66 8.66 0.32 0.31 8.15 9.19 1.00 1136 1244
## m1[24] 7.98 7.98 0.28 0.27 7.53 8.45 1.00 1418 1074
## m1[25] 0.89 0.90 0.34 0.33 0.30 1.40 1.01 386 466
## m1[26] 7.73 7.73 0.27 0.26 7.30 8.18 1.00 1571 1003
## m1[27] 3.23 3.24 0.22 0.22 2.87 3.58 1.01 498 657
## m1[28] 3.39 3.39 0.22 0.21 3.04 3.73 1.01 522 665
## m1[29] 8.27 8.27 0.30 0.29 7.80 8.77 1.00 1278 1179
## m1[30] 2.88 2.88 0.23 0.23 2.48 3.25 1.01 448 591
y1=mcmc$draws('y1')
smy=summary(y1)
tb$y1=smy$median
tb$l5=smy$q5
tb$u5=smy$q95
par(mfrow=c(1,1))
plot(tb$y,tb$y1,xlab='obs.',ylab='prd.')
abline(0,1)
qplot(data=tb,y,y1,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1)
par(mfrow=c(1,2))
hist(tb$y-tb$y1,xlab='obs.-prd.',main='residual')
density(tb$y-tb$y1) |> plot(xlab='obs.-prd.',main='residual')
tb=arrange(tb,y1)
par(mfrow=c(1,1))
ylim=c(min(tb$l5),max(tb$u5))
plot(tb$y,ylim=ylim,ylab='y',col='red')
par(new=T)
plot(tb$y1,ylim=ylim,ylab='',col='black')
par(new=T)
plot(tb$l5,ylim=ylim,ylab='',type='l',col='gray')
par(new=T)
plot(tb$u5,ylim=ylim,ylab='',type='l',col='gray')
qplot(1:n,y,data=tb,col=I('red'))+
geom_point(aes(x=1:n,y=y1),tb,col='black')+
geom_line(aes(x=1:n,y=l5),tb,col='gray')+
geom_line(aes(x=1:n,y=u5),tb,col='gray')
multiple regression
tb=tibble(x1=runif(n,0,9),x2=runif(n,0,9),
y=rnorm(n,x1-x2,1))
f=formula(y~x1+x2)
par(mfrow=c(1,2))
plot(tb$x1,tb$y)
plot(tb$x2,tb$y)
X=model.matrix(f,tb)
data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X)
mle=mdl$optimize(data=data)## Initial log joint probability = -4204.6
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 25 -14.3992 5.25865e-05 0.000541205 1 1 54
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -14.40
## b[1] 0.39
## b[2] 0.93
## b[3] -1.05
## s 0.98
## y1[1] 3.42
## y1[2] 6.67
## y1[3] -3.01
## y1[4] 6.78
## y1[5] -0.42
## y1[6] 3.68
## y1[7] -4.01
## y1[8] 1.52
## y1[9] 4.15
## y1[10] -5.76
## y1[11] -1.92
## y1[12] -2.42
## y1[13] 1.00
## y1[14] 5.72
## y1[15] -4.48
## y1[16] 3.54
## y1[17] 1.75
## y1[18] -2.47
## y1[19] 5.21
## y1[20] -2.72
## y1[21] 0.27
## y1[22] 1.25
## y1[23] -2.80
## y1[24] 4.40
## y1[25] -3.94
## y1[26] 0.53
## y1[27] -0.15
## y1[28] -4.72
## y1[29] 3.05
## y1[30] -6.67
## m1[1] 3.23
## m1[2] 5.44
## m1[3] -2.45
## m1[4] 7.59
## m1[5] -1.23
## m1[6] 4.17
## m1[7] -5.68
## m1[8] 1.46
## m1[9] 3.94
## m1[10] -5.36
## m1[11] -1.49
## m1[12] -2.82
## m1[13] 1.84
## m1[14] 5.50
## m1[15] -3.57
## m1[16] 3.31
## m1[17] 1.24
## m1[18] -2.68
## m1[19] 5.61
## m1[20] -2.56
## m1[21] 0.13
## m1[22] -0.18
## m1[23] -2.46
## m1[24] 3.92
## m1[25] -3.29
## m1[26] -0.32
## m1[27] -0.87
## m1[28] -3.23
## m1[29] 4.06
## m1[30] -6.36mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc,exc='m',ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -16.62 -16.29 1.56 1.38 -19.42 -14.80 1.00 563 791
## b[1] 0.33 0.34 0.53 0.50 -0.58 1.17 1.01 348 284
## b[2] 0.94 0.94 0.08 0.08 0.81 1.07 1.01 652 730
## b[3] -1.04 -1.04 0.07 0.07 -1.16 -0.93 1.00 799 767
## s 1.09 1.07 0.16 0.15 0.86 1.38 1.00 983 958
## y1[1] 3.24 3.27 1.13 1.09 1.29 5.06 1.00 1845 1992
## y1[2] 5.45 5.44 1.11 1.09 3.61 7.30 1.00 1989 1677
## y1[3] -2.53 -2.51 1.20 1.20 -4.45 -0.58 1.00 1532 1712
## y1[4] 7.56 7.57 1.21 1.12 5.59 9.60 1.00 1884 1854
## y1[5] -1.25 -1.22 1.14 1.10 -3.10 0.59 1.00 1575 1842
## y1[6] 4.20 4.25 1.13 1.11 2.32 6.07 1.00 1621 1689
## y1[7] -5.71 -5.67 1.15 1.14 -7.60 -3.91 1.00 1790 1931
## y1[8] 1.40 1.42 1.17 1.13 -0.52 3.30 1.00 1275 1524
## y1[9] 3.91 3.88 1.18 1.11 2.08 5.83 1.00 1935 1856
## y1[10] -5.35 -5.37 1.16 1.11 -7.25 -3.43 1.00 2018 1812
## y1[11] -1.49 -1.50 1.12 1.10 -3.33 0.30 1.00 1901 1697
## y1[12] -2.82 -2.82 1.13 1.09 -4.70 -1.02 1.00 1630 1895
## y1[13] 1.83 1.87 1.18 1.10 -0.22 3.70 1.00 1880 1953
## y1[14] 5.47 5.44 1.18 1.14 3.63 7.45 1.00 1916 1973
## y1[15] -3.58 -3.60 1.21 1.20 -5.54 -1.57 1.00 1753 1974
## y1[16] 3.31 3.31 1.15 1.12 1.40 5.20 1.00 1551 1572
## y1[17] 1.27 1.25 1.10 1.07 -0.52 3.08 1.00 2020 2000
## y1[18] -2.65 -2.68 1.15 1.13 -4.53 -0.82 1.00 1813 1862
## y1[19] 5.66 5.70 1.17 1.12 3.70 7.60 1.00 2002 1922
## y1[20] -2.52 -2.50 1.14 1.11 -4.35 -0.67 1.00 1988 1880
## y1[21] 0.09 0.09 1.16 1.10 -1.83 1.94 1.00 1719 1935
## y1[22] -0.16 -0.19 1.20 1.20 -2.07 1.79 1.00 1740 1771
## y1[23] -2.45 -2.45 1.17 1.19 -4.32 -0.57 1.00 2083 1968
## y1[24] 3.88 3.89 1.17 1.13 1.98 5.74 1.00 1781 1881
## y1[25] -3.34 -3.31 1.13 1.11 -5.20 -1.51 1.00 1829 1948
## y1[26] -0.35 -0.33 1.10 1.11 -2.10 1.39 1.00 1839 1404
## y1[27] -0.85 -0.87 1.17 1.16 -2.68 1.05 1.00 1771 1894
## y1[28] -3.27 -3.27 1.14 1.10 -5.16 -1.42 1.00 2029 1893
## y1[29] 4.06 4.10 1.12 1.11 2.22 5.90 1.00 1883 1926
## y1[30] -6.39 -6.39 1.18 1.15 -8.32 -4.48 1.00 1743 1736
## m1[1] 3.22 3.22 0.29 0.27 2.74 3.69 1.00 629 739
## m1[2] 5.44 5.45 0.36 0.35 4.85 6.01 1.00 1600 1019
## m1[3] -2.49 -2.49 0.40 0.39 -3.18 -1.87 1.01 386 380
## m1[4] 7.60 7.60 0.47 0.45 6.81 8.34 1.00 1841 1278
## m1[5] -1.27 -1.26 0.34 0.34 -1.84 -0.74 1.01 372 382
## m1[6] 4.15 4.15 0.35 0.34 3.56 4.70 1.00 566 604
## m1[7] -5.69 -5.69 0.38 0.36 -6.32 -5.06 1.00 1609 1392
## m1[8] 1.42 1.42 0.39 0.36 0.73 2.04 1.01 359 343
## m1[9] 3.94 3.95 0.30 0.29 3.44 4.43 1.00 1735 1348
## m1[10] -5.36 -5.36 0.36 0.35 -5.97 -4.75 1.00 2130 1380
## m1[11] -1.49 -1.49 0.24 0.23 -1.91 -1.08 1.00 2153 1440
## m1[12] -2.82 -2.82 0.30 0.30 -3.31 -2.32 1.00 1743 1714
## m1[13] 1.86 1.87 0.34 0.34 1.31 2.39 1.00 1189 1389
## m1[14] 5.50 5.51 0.37 0.35 4.90 6.08 1.00 1697 1088
## m1[15] -3.60 -3.60 0.36 0.36 -4.22 -3.02 1.01 489 587
## m1[16] 3.28 3.28 0.39 0.38 2.61 3.90 1.01 409 344
## m1[17] 1.25 1.25 0.23 0.23 0.86 1.62 1.00 1957 1671
## m1[18] -2.68 -2.68 0.29 0.28 -3.15 -2.20 1.00 1936 1647
## m1[19] 5.61 5.61 0.37 0.36 5.00 6.18 1.00 1595 1052
## m1[20] -2.54 -2.54 0.38 0.37 -3.16 -1.91 1.00 1047 1245
## m1[21] 0.09 0.10 0.37 0.36 -0.57 0.68 1.01 354 327
## m1[22] -0.15 -0.15 0.46 0.45 -0.90 0.57 1.01 728 915
## m1[23] -2.44 -2.45 0.34 0.33 -3.00 -1.88 1.00 1267 1451
## m1[24] 3.90 3.90 0.32 0.30 3.38 4.41 1.00 661 825
## m1[25] -3.31 -3.30 0.29 0.28 -3.79 -2.84 1.00 873 1190
## m1[26] -0.34 -0.33 0.24 0.24 -0.73 0.02 1.00 492 555
## m1[27] -0.84 -0.84 0.42 0.42 -1.53 -0.17 1.01 801 1081
## m1[28] -3.25 -3.24 0.27 0.26 -3.71 -2.80 1.00 1419 1515
## m1[29] 4.07 4.08 0.32 0.32 3.53 4.59 1.00 1879 1735
## m1[30] -6.39 -6.39 0.43 0.42 -7.10 -5.67 1.00 996 1238
y1=mcmc$draws('y1')
smy=summary(y1)
tb$y1=smy$median
tb$l5=smy$q5
tb$u5=smy$q95
par(mfrow=c(1,1))
plot(tb$y,tb$y1,xlab='obs.',ylab='prd.')
abline(0,1)
qplot(data=tb,y,y1,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1)
par(mfrow=c(1,2))
hist(tb$y-tb$y1,xlab='obs.-prd.',main='residual')
density(tb$y-tb$y1) |> plot(xlab='obs.-prd.',main='residual')
tb=arrange(tb,y1)
par(mfrow=c(1,1))
ylim=c(min(tb$l5),max(tb$u5))
plot(tb$y,ylim=ylim,ylab='y',col='red')
par(new=T)
plot(tb$y1,ylim=ylim,ylab='',col='black')
par(new=T)
plot(tb$l5,ylim=ylim,ylab='',type='l',col='gray')
par(new=T)
plot(tb$u5,ylim=ylim,ylab='',type='l',col='gray')
qplot(1:n,y,data=tb,col=I('red'))+
geom_point(aes(x=1:n,y=y1),tb,col='black')+
geom_line(aes(x=1:n,y=l5),tb,col='gray')+
geom_line(aes(x=1:n,y=u5),tb,col='gray')
ANOVA
tb=tibble(c=sample(c('a','b','c'),n,replace=T),
y=rnorm(n,(c=='b')*2-(c=='c')*2,1))
f=formula(y~c)
par(mfrow=c(1,1))
boxplot(y~c,tb)
qplot(data=tb,c,y,geom='boxplot')
X=model.matrix(f,tb)
data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X)
mle=mdl$optimize(data=data)## Initial log joint probability = -29.1885
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 18 -12.4727 5.90482e-05 6.64483e-05 0.8902 0.8902 21
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -12.47
## b[1] 0.20
## b[2] 1.72
## b[3] -2.18
## s 0.92
## y1[1] -2.72
## y1[2] 0.51
## y1[3] -0.40
## y1[4] 0.74
## y1[5] 0.36
## y1[6] -0.08
## y1[7] 1.91
## y1[8] 2.80
## y1[9] 2.74
## y1[10] 0.24
## y1[11] 0.06
## y1[12] 3.02
## y1[13] -0.92
## y1[14] 1.17
## y1[15] -2.18
## y1[16] 1.40
## y1[17] 0.09
## y1[18] 2.94
## y1[19] -0.20
## y1[20] 0.74
## y1[21] 1.68
## y1[22] -0.47
## y1[23] 0.26
## y1[24] 0.39
## y1[25] 3.70
## y1[26] -1.97
## y1[27] -1.66
## y1[28] 0.91
## y1[29] 0.47
## y1[30] 1.82
## m1[1] -1.98
## m1[2] 0.20
## m1[3] 0.20
## m1[4] 0.20
## m1[5] 0.20
## m1[6] 0.20
## m1[7] 1.92
## m1[8] 1.92
## m1[9] 1.92
## m1[10] 0.20
## m1[11] 0.20
## m1[12] 1.92
## m1[13] 0.20
## m1[14] 1.92
## m1[15] -1.98
## m1[16] 1.92
## m1[17] 0.20
## m1[18] 1.92
## m1[19] -1.98
## m1[20] 1.92
## m1[21] 1.92
## m1[22] 0.20
## m1[23] 0.20
## m1[24] -1.98
## m1[25] 1.92
## m1[26] -1.98
## m1[27] -1.98
## m1[28] 1.92
## m1[29] 0.20
## m1[30] 1.92mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc,exc='m',ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -14.71 -14.34 1.53 1.28 -17.67 -12.90 1.00 785 1123
## b[1] 0.19 0.20 0.30 0.31 -0.32 0.67 1.00 867 1175
## b[2] 1.73 1.71 0.42 0.40 1.08 2.48 1.00 883 768
## b[3] -2.18 -2.18 0.50 0.48 -3.02 -1.38 1.00 682 892
## s 1.01 1.00 0.15 0.14 0.80 1.27 1.00 2321 1565
## y1[1] -1.96 -1.96 1.12 1.10 -3.74 -0.02 1.00 1959 1871
## y1[2] 0.16 0.15 1.09 1.05 -1.62 1.93 1.00 1816 1751
## y1[3] 0.20 0.22 1.08 1.06 -1.64 1.99 1.00 1815 2000
## y1[4] 0.18 0.20 1.06 1.06 -1.54 1.93 1.00 1862 1753
## y1[5] 0.21 0.24 1.07 1.02 -1.61 1.94 1.00 1777 1638
## y1[6] 0.20 0.19 1.05 1.00 -1.47 1.95 1.01 1510 1712
## y1[7] 1.89 1.85 1.12 1.06 0.09 3.73 1.00 1936 1996
## y1[8] 1.91 1.90 1.06 1.07 0.18 3.64 1.00 1994 1980
## y1[9] 1.96 1.97 1.06 1.05 0.21 3.69 1.00 1877 1764
## y1[10] 0.18 0.19 1.04 0.98 -1.60 1.98 1.00 1654 1785
## y1[11] 0.20 0.20 1.05 1.01 -1.49 1.92 1.00 1598 1955
## y1[12] 1.90 1.88 1.08 1.03 0.08 3.67 1.00 2085 1875
## y1[13] 0.20 0.20 1.07 1.08 -1.58 1.90 1.00 2031 1921
## y1[14] 1.87 1.90 1.10 1.03 0.08 3.60 1.00 1994 1805
## y1[15] -1.98 -1.99 1.08 1.05 -3.82 -0.22 1.00 1754 1615
## y1[16] 1.92 1.90 1.08 1.05 0.17 3.68 1.00 2115 1554
## y1[17] 0.18 0.19 1.09 1.04 -1.61 2.00 1.00 1484 1807
## y1[18] 1.95 1.97 1.09 1.02 0.10 3.78 1.00 1755 1929
## y1[19] -2.02 -2.02 1.14 1.10 -3.92 -0.15 1.00 1651 1958
## y1[20] 1.92 1.91 1.08 1.03 0.16 3.71 1.00 1610 1795
## y1[21] 1.91 1.92 1.08 1.02 0.13 3.64 1.00 2075 1957
## y1[22] 0.20 0.22 1.07 1.06 -1.56 1.94 1.00 1790 1804
## y1[23] 0.20 0.20 1.07 1.06 -1.56 1.98 1.00 1672 1643
## y1[24] -2.00 -2.01 1.12 1.11 -3.83 -0.18 1.00 1915 1772
## y1[25] 1.94 1.94 1.05 1.02 0.17 3.60 1.00 1681 1908
## y1[26] -1.97 -1.94 1.11 1.12 -3.83 -0.14 1.00 1677 1881
## y1[27] -1.99 -1.98 1.07 1.07 -3.75 -0.30 1.00 1818 1690
## y1[28] 1.93 1.93 1.04 1.01 0.19 3.68 1.00 1970 1953
## y1[29] 0.21 0.21 1.06 1.07 -1.54 1.94 1.00 1516 1456
## y1[30] 1.96 1.93 1.08 1.06 0.18 3.75 1.00 1753 1871
## m1[1] -1.98 -1.98 0.41 0.39 -2.65 -1.32 1.00 1092 1331
## m1[2] 0.19 0.20 0.30 0.31 -0.32 0.67 1.00 867 1175
## m1[3] 0.19 0.20 0.30 0.31 -0.32 0.67 1.00 867 1175
## m1[4] 0.19 0.20 0.30 0.31 -0.32 0.67 1.00 867 1175
## m1[5] 0.19 0.20 0.30 0.31 -0.32 0.67 1.00 867 1175
## m1[6] 0.19 0.20 0.30 0.31 -0.32 0.67 1.00 867 1175
## m1[7] 1.92 1.92 0.30 0.28 1.46 2.41 1.00 1508 1272
## m1[8] 1.92 1.92 0.30 0.28 1.46 2.41 1.00 1508 1272
## m1[9] 1.92 1.92 0.30 0.28 1.46 2.41 1.00 1508 1272
## m1[10] 0.19 0.20 0.30 0.31 -0.32 0.67 1.00 867 1175
## m1[11] 0.19 0.20 0.30 0.31 -0.32 0.67 1.00 867 1175
## m1[12] 1.92 1.92 0.30 0.28 1.46 2.41 1.00 1508 1272
## m1[13] 0.19 0.20 0.30 0.31 -0.32 0.67 1.00 867 1175
## m1[14] 1.92 1.92 0.30 0.28 1.46 2.41 1.00 1508 1272
## m1[15] -1.98 -1.98 0.41 0.39 -2.65 -1.32 1.00 1092 1331
## m1[16] 1.92 1.92 0.30 0.28 1.46 2.41 1.00 1508 1272
## m1[17] 0.19 0.20 0.30 0.31 -0.32 0.67 1.00 867 1175
## m1[18] 1.92 1.92 0.30 0.28 1.46 2.41 1.00 1508 1272
## m1[19] -1.98 -1.98 0.41 0.39 -2.65 -1.32 1.00 1092 1331
## m1[20] 1.92 1.92 0.30 0.28 1.46 2.41 1.00 1508 1272
## m1[21] 1.92 1.92 0.30 0.28 1.46 2.41 1.00 1508 1272
## m1[22] 0.19 0.20 0.30 0.31 -0.32 0.67 1.00 867 1175
## m1[23] 0.19 0.20 0.30 0.31 -0.32 0.67 1.00 867 1175
## m1[24] -1.98 -1.98 0.41 0.39 -2.65 -1.32 1.00 1092 1331
## m1[25] 1.92 1.92 0.30 0.28 1.46 2.41 1.00 1508 1272
## m1[26] -1.98 -1.98 0.41 0.39 -2.65 -1.32 1.00 1092 1331
## m1[27] -1.98 -1.98 0.41 0.39 -2.65 -1.32 1.00 1092 1331
## m1[28] 1.92 1.92 0.30 0.28 1.46 2.41 1.00 1508 1272
## m1[29] 0.19 0.20 0.30 0.31 -0.32 0.67 1.00 867 1175
## m1[30] 1.92 1.92 0.30 0.28 1.46 2.41 1.00 1508 1272
y1=mcmc$draws('y1')
smy=summary(y1)
tb$y1=smy$median
tb$l5=smy$q5
tb$u5=smy$q95
lv=c(0,1,2)
par(mfrow=c(1,1))
plot(tb$y,tb$y1,pch=lv[factor(tb$c)],xlab='obs.',ylab='prd.')
qplot(data=tb,y,y1,shape=c,size=I(2),xlab='obs.',ylab='prd.')
plot(tb$y,tb$y1,col=1+lv[factor(tb$c)],xlab='obs.',ylab='prd.')
qplot(data=tb,y,y1,col=c,xlab='obs.',ylab='prd.')
par(mfrow=c(1,2))
hist(tb$y-tb$y1,xlab='obs.-prd.',main='residual')
density(tb$y-tb$y1) |> plot(xlab='obs.-prd.',main='residual')
tb=arrange(tb,y1)
par(mfrow=c(1,1))
ylim=c(min(tb$l5),max(tb$u5))
plot(tb$y,ylim=ylim,ylab='y',col='red',pch=lv[factor(tb$c)])
par(new=T)
plot(tb$y1,ylim=ylim,ylab='',col='black',pch=lv[factor(tb$c)])
par(new=T)
plot(tb$l5,ylim=ylim,ylab='',type='l',col='gray')
par(new=T)
plot(tb$u5,ylim=ylim,ylab='',type='l',col='gray')
qplot(data=tb,1:n,y,col=c)+
geom_point(aes(x=1:n,y=y1),tb,col='black')+
geom_line(aes(x=1:n,y=l5),tb,col='gray')+
geom_line(aes(x=1:n,y=u5),tb,col='gray')
ANCOVA
tb=tibble(x=runif(n,0,9),c=sample(c('a','b','c'),n,replace=T),
y=rnorm(n,x+(c=='b')*2-(c=='c')*2,1))
f=formula(y~x+c)
par(mfrow=c(1,1))
#plot(tb$x1,tb$y,pch=tb$c1)
lv=c(0,1,2)
plot(tb$x,tb$y,pch=lv[factor(tb$c)])
qplot(data=tb,x,y,shape=c,size=I(2))
plot(tb$x,tb$y,col=1+lv[factor(tb$c)])
qplot(data=tb,x,y,col=c)
X=model.matrix(f,tb)
data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X)
mle=mdl$optimize(data=data)## Initial log joint probability = -152.843
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 40 -8.25775 5.84488e-05 0.00125024 1 1 49
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -8.26
## b[1] -0.08
## b[2] 0.95
## b[3] 2.86
## b[4] -1.84
## s 0.80
## y1[1] 4.76
## y1[2] -1.93
## y1[3] 0.72
## y1[4] 0.93
## y1[5] 5.97
## y1[6] 8.09
## y1[7] 10.46
## y1[8] 7.13
## y1[9] 5.33
## y1[10] 6.06
## y1[11] -0.12
## y1[12] -1.08
## y1[13] 6.56
## y1[14] 5.25
## y1[15] 2.12
## y1[16] 6.44
## y1[17] 4.81
## y1[18] 0.49
## y1[19] 5.53
## y1[20] 4.02
## y1[21] 8.40
## y1[22] 0.32
## y1[23] 4.53
## y1[24] 4.66
## y1[25] 4.79
## y1[26] 4.16
## y1[27] 7.74
## y1[28] 8.27
## y1[29] 5.39
## y1[30] 2.88
## m1[1] 4.30
## m1[2] -1.68
## m1[3] 1.11
## m1[4] 0.57
## m1[5] 4.59
## m1[6] 8.26
## m1[7] 11.17
## m1[8] 8.04
## m1[9] 5.40
## m1[10] 6.54
## m1[11] -0.12
## m1[12] -0.64
## m1[13] 6.19
## m1[14] 5.64
## m1[15] 1.70
## m1[16] 6.49
## m1[17] 4.00
## m1[18] 1.01
## m1[19] 4.75
## m1[20] 2.92
## m1[21] 7.76
## m1[22] 1.05
## m1[23] 4.55
## m1[24] 4.51
## m1[25] 4.23
## m1[26] 3.76
## m1[27] 7.58
## m1[28] 8.59
## m1[29] 5.39
## m1[30] 3.28mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc,exc='m',ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -11.19 -10.84 1.72 1.55 -14.42 -9.09 1.01 751 1048
## b[1] -0.09 -0.09 0.34 0.33 -0.64 0.49 1.00 925 949
## b[2] 0.95 0.95 0.06 0.06 0.85 1.04 1.00 1606 1661
## b[3] 2.86 2.85 0.37 0.35 2.28 3.47 1.01 814 1121
## b[4] -1.81 -1.81 0.45 0.44 -2.56 -1.07 1.01 799 857
## s 0.91 0.89 0.13 0.12 0.71 1.14 1.00 2143 1743
## y1[1] 4.30 4.30 0.94 0.93 2.72 5.79 1.00 1959 1990
## y1[2] -1.66 -1.68 1.05 1.02 -3.35 0.09 1.00 1692 1814
## y1[3] 1.09 1.11 0.98 0.92 -0.56 2.73 1.00 1511 1868
## y1[4] 0.53 0.54 0.96 0.91 -1.04 2.06 1.00 1938 1985
## y1[5] 4.61 4.62 0.96 0.92 3.00 6.22 1.00 1837 1903
## y1[6] 8.24 8.23 0.96 0.93 6.74 9.83 1.00 1770 2031
## y1[7] 11.18 11.19 1.00 1.00 9.58 12.80 1.00 1951 1917
## y1[8] 8.01 8.03 0.96 0.92 6.43 9.54 1.00 1850 1894
## y1[9] 5.35 5.35 0.98 0.97 3.73 6.92 1.00 1865 1902
## y1[10] 6.53 6.52 1.01 1.00 4.90 8.14 1.00 1933 1958
## y1[11] -0.11 -0.10 1.01 0.96 -1.75 1.50 1.00 1501 1866
## y1[12] -0.64 -0.64 1.00 0.97 -2.25 0.98 1.00 1797 1855
## y1[13] 6.19 6.17 0.94 0.90 4.69 7.77 1.00 1642 1728
## y1[14] 5.67 5.66 1.02 0.98 3.95 7.32 1.00 2042 1954
## y1[15] 1.70 1.70 0.95 0.91 0.14 3.23 1.00 1759 1893
## y1[16] 6.46 6.46 0.97 0.95 4.86 8.01 1.00 2145 1963
## y1[17] 4.00 3.99 0.94 0.89 2.47 5.52 1.00 1891 1842
## y1[18] 0.98 0.95 1.00 1.02 -0.63 2.63 1.00 1733 1584
## y1[19] 4.77 4.80 0.95 0.89 3.24 6.30 1.00 1939 1807
## y1[20] 2.93 2.90 0.98 0.93 1.37 4.59 1.00 1874 1855
## y1[21] 7.72 7.72 0.95 0.93 6.13 9.32 1.00 1923 1969
## y1[22] 1.04 1.04 0.93 0.92 -0.51 2.56 1.00 1846 1906
## y1[23] 4.52 4.55 0.98 0.91 2.92 6.14 1.00 2052 1919
## y1[24] 4.51 4.53 0.97 0.95 2.94 6.10 1.00 1733 1863
## y1[25] 4.19 4.20 0.94 0.94 2.66 5.71 1.00 1929 1930
## y1[26] 3.73 3.71 0.93 0.93 2.18 5.26 1.00 1985 1933
## y1[27] 7.56 7.57 0.98 0.94 5.97 9.21 1.00 1809 1848
## y1[28] 8.60 8.61 0.99 0.97 6.99 10.18 1.00 1960 1920
## y1[29] 5.36 5.37 0.94 0.94 3.78 6.86 1.00 1909 1940
## y1[30] 3.32 3.34 0.99 0.97 1.63 4.94 1.00 1803 1766
## m1[1] 4.29 4.28 0.30 0.29 3.80 4.80 1.00 1471 1122
## m1[2] -1.67 -1.67 0.46 0.43 -2.42 -0.89 1.01 1077 1000
## m1[3] 1.09 1.09 0.30 0.29 0.63 1.59 1.00 875 1144
## m1[4] 0.55 0.55 0.32 0.31 0.05 1.08 1.00 891 1031
## m1[5] 4.58 4.58 0.25 0.25 4.16 4.98 1.01 992 1322
## m1[6] 8.25 8.25 0.31 0.29 7.74 8.75 1.00 1228 1410
## m1[7] 11.15 11.15 0.42 0.39 10.46 11.84 1.00 1303 1267
## m1[8] 8.03 8.02 0.30 0.28 7.54 8.53 1.00 1227 1411
## m1[9] 5.38 5.38 0.28 0.27 4.94 5.85 1.00 1372 1288
## m1[10] 6.55 6.56 0.46 0.45 5.78 7.27 1.01 1405 1407
## m1[11] -0.11 -0.11 0.42 0.38 -0.80 0.58 1.01 1056 885
## m1[12] -0.63 -0.63 0.43 0.40 -1.34 0.09 1.01 1060 826
## m1[13] 6.17 6.17 0.28 0.26 5.73 6.63 1.00 1265 1306
## m1[14] 5.66 5.67 0.43 0.42 4.93 6.32 1.01 1332 1292
## m1[15] 1.69 1.69 0.28 0.26 1.25 2.15 1.00 852 1034
## m1[16] 6.48 6.49 0.30 0.30 5.96 6.96 1.01 1287 1254
## m1[17] 3.99 4.00 0.25 0.24 3.58 4.38 1.01 926 1203
## m1[18] 0.99 0.99 0.30 0.29 0.52 1.50 1.00 881 1114
## m1[19] 4.74 4.75 0.26 0.25 4.32 5.15 1.01 1017 1271
## m1[20] 2.90 2.89 0.34 0.33 2.33 3.49 1.00 1570 1211
## m1[21] 7.75 7.76 0.36 0.35 7.15 8.32 1.01 1503 1321
## m1[22] 1.03 1.03 0.30 0.29 0.56 1.53 1.00 878 1117
## m1[23] 4.53 4.52 0.29 0.28 4.06 5.03 1.00 1452 1192
## m1[24] 4.50 4.49 0.30 0.28 4.02 5.00 1.00 1454 1192
## m1[25] 4.22 4.22 0.25 0.24 3.80 4.61 1.01 949 1251
## m1[26] 3.75 3.75 0.25 0.24 3.34 4.14 1.01 904 1225
## m1[27] 7.57 7.59 0.35 0.34 6.99 8.13 1.01 1464 1326
## m1[28] 8.58 8.58 0.31 0.30 8.05 9.10 1.00 1238 1331
## m1[29] 5.38 5.39 0.27 0.26 4.92 5.81 1.01 1109 1282
## m1[30] 3.29 3.31 0.39 0.37 2.65 3.91 1.01 1150 994
y1=mcmc$draws('y1')
smy=summary(y1)
tb$y1=smy$median
tb$l5=smy$q5
tb$u5=smy$q95
lv=c(0,1,2)
par(mfrow=c(1,1))
plot(tb$y,tb$y1,pch=lv[factor(tb$c)],xlab='obs.',ylab='prd.')
abline(0,1)
plot(tb$y,tb$y1,col=1+lv[factor(tb$c)],xlab='obs.',ylab='prd.')
abline(0,1)
qplot(data=tb,y,y1,col=c,xlab='obs.',ylab='prd.')+
geom_abline(intercept=0,slope=1)
par(mfrow=c(1,2))
hist(tb$y-tb$y1,xlab='obs.-prd.',main='residual')
density(tb$y-tb$y1) |> plot(xlab='obs.-prd.',main='residual')
tb=arrange(tb,y1)
par(mfrow=c(1,1))
ylim=c(min(tb$l5),max(tb$u5))
plot(tb$y,ylim=ylim,ylab='y',col='red',pch=lv[factor(tb$c)])
par(new=T)
plot(tb$y1,ylim=ylim,ylab='',col='black',pch=lv[factor(tb$c)])
par(new=T)
plot(tb$l5,ylim=ylim,ylab='',type='l',col='gray')
par(new=T)
plot(tb$u5,ylim=ylim,ylab='',type='l',col='gray')
qplot(data=tb,1:n,y,col=c)+
geom_point(aes(x=1:n,y=y1),tb,col='black')+
geom_line(aes(x=1:n,y=l5),tb,col='gray')+
geom_line(aes(x=1:n,y=u5),tb,col='gray')
ex6
generalized linear regression
poisson regression
objective variable [0,Infinity), integer. variance of error is near to mean
(normal linear regression, correlation is near to 1,-1, variance of error is constant)
# of samples is N,
log li=sum(bj*xji),j=[0,K],i=[1,N]
yi~Po(li)
or
li=sum(bj xji),j=[0,k]
yi~Po(exp li)
when xj -> xj +1, y -> y* exp bj ex6-1.stan
data{
int N;
int K;
array[N] int y;
matrix[N,K] X;
}
parameters{
vector[K] b;
}
model{
vector[N] l=X*b;
y~poisson_log(l);
}
generated quantities{
array[N] int y1;
vector[N] l1=X*b;
for(i in 1:N){
y1[i]=poisson_log_rng(l1[i]);
}
}n=30
tb=tibble(x=runif(n,-1,1),c=sample(c('a','b'),n,replace=T),
y=rpois(n,exp(x+(c=='b')*0.5)))
f=formula(y~x+c)
qplot(data=tb,x,y,col=c)
glm(f,tb,family='poisson')##
## Call: glm(formula = f, family = "poisson", data = tb)
##
## Coefficients:
## (Intercept) x cb
## 0.0658 1.0269 0.7378
##
## Degrees of Freedom: 29 Total (i.e. Null); 27 Residual
## Null Deviance: 48.5
## Residual Deviance: 26.7 AIC: 90.2X=model.matrix(f,tb)
data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X)
mdl=cmdstan_model('./ex6-1.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -61.439
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 8 -13.0485 0.000132489 0.000119888 1 1 11
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -13.05
## b[1] 0.07
## b[2] 1.03
## b[3] 0.74
## y1[1] 0.00
## y1[2] 0.00
## y1[3] 2.00
## y1[4] 5.00
## y1[5] 0.00
## y1[6] 2.00
## y1[7] 0.00
## y1[8] 7.00
## y1[9] 2.00
## y1[10] 2.00
## y1[11] 1.00
## y1[12] 1.00
## y1[13] 0.00
## y1[14] 1.00
## y1[15] 1.00
## y1[16] 0.00
## y1[17] 0.00
## y1[18] 1.00
## y1[19] 2.00
## y1[20] 4.00
## y1[21] 3.00
## y1[22] 1.00
## y1[23] 1.00
## y1[24] 0.00
## y1[25] 6.00
## y1[26] 1.00
## y1[27] 3.00
## y1[28] 2.00
## y1[29] 7.00
## y1[30] 1.00
## l1[1] -0.09
## l1[2] -0.17
## l1[3] 0.27
## l1[4] 1.66
## l1[5] -0.23
## l1[6] 0.98
## l1[7] 0.62
## l1[8] 1.04
## l1[9] 0.89
## l1[10] 0.10
## l1[11] 0.16
## l1[12] 0.03
## l1[13] -0.64
## l1[14] 0.07
## l1[15] 0.23
## l1[16] -0.14
## l1[17] -0.18
## l1[18] 0.53
## l1[19] 0.33
## l1[20] 1.63
## l1[21] 0.04
## l1[22] -0.90
## l1[23] 0.23
## l1[24] -0.77
## l1[25] 0.85
## l1[26] 0.55
## l1[27] 0.60
## l1[28] 0.64
## l1[29] 1.41
## l1[30] -0.06mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc,exc='l1',ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -14.53 -14.30 1.11 0.98 -16.54 -13.25 1.00 1024 1083
## b[1] 0.04 0.04 0.21 0.22 -0.30 0.38 1.00 1110 1283
## b[2] 1.03 1.02 0.25 0.26 0.63 1.44 1.00 1289 1447
## b[3] 0.74 0.74 0.28 0.29 0.27 1.18 1.00 1073 1312
## y1[1] 0.96 1.00 1.05 1.48 0.00 3.00 1.00 1755 1590
## y1[2] 0.82 1.00 0.93 1.48 0.00 3.00 1.00 2041 1929
## y1[3] 1.34 1.00 1.15 1.48 0.00 3.00 1.00 1937 1960
## y1[4] 5.20 5.00 2.63 2.97 1.00 10.00 1.00 1732 1962
## y1[5] 0.81 1.00 0.92 1.48 0.00 3.00 1.00 1888 1862
## y1[6] 2.67 2.00 1.71 1.48 0.00 6.00 1.00 1933 2055
## y1[7] 1.78 2.00 1.40 1.48 0.00 4.00 1.00 2060 2048
## y1[8] 2.84 3.00 1.76 1.48 0.00 6.00 1.00 1923 1968
## y1[9] 2.33 2.00 1.61 1.48 0.00 5.00 1.00 1891 1881
## y1[10] 1.11 1.00 1.09 1.48 0.00 3.00 1.00 2036 1987
## y1[11] 1.20 1.00 1.15 1.48 0.00 3.00 1.00 1925 1893
## y1[12] 1.01 1.00 1.04 1.48 0.00 3.00 1.00 1812 1577
## y1[13] 0.57 0.00 0.80 0.00 0.00 2.00 1.00 1842 2025
## y1[14] 1.12 1.00 1.08 1.48 0.00 3.00 1.00 1863 1843
## y1[15] 1.23 1.00 1.12 1.48 0.00 3.00 1.00 1884 1909
## y1[16] 0.86 1.00 0.95 1.48 0.00 3.00 1.00 2002 1982
## y1[17] 0.90 1.00 0.99 1.48 0.00 3.00 1.00 2004 1813
## y1[18] 1.69 1.00 1.34 1.48 0.00 4.00 1.00 2030 1915
## y1[19] 1.36 1.00 1.24 1.48 0.00 4.00 1.00 2005 1990
## y1[20] 5.10 5.00 2.56 2.97 1.00 10.00 1.00 1902 2084
## y1[21] 1.00 1.00 1.00 1.48 0.00 3.00 1.00 2086 2021
## y1[22] 0.40 0.00 0.66 0.00 0.00 2.00 1.00 1791 1817
## y1[23] 1.23 1.00 1.11 1.48 0.00 3.00 1.00 1890 1916
## y1[24] 0.45 0.00 0.71 0.00 0.00 2.00 1.00 1905 1948
## y1[25] 2.29 2.00 1.61 1.48 0.00 5.00 1.00 2111 1993
## y1[26] 1.71 2.00 1.38 1.48 0.00 4.00 1.00 2126 1936
## y1[27] 1.88 2.00 1.46 1.48 0.00 5.00 1.00 1892 1949
## y1[28] 1.86 2.00 1.41 1.48 0.00 4.00 1.00 1887 1984
## y1[29] 4.19 4.00 2.30 2.97 1.00 8.00 1.00 1616 1904
## y1[30] 0.90 1.00 0.97 1.48 0.00 3.00 1.00 2042 1935
## l1[1] -0.11 -0.10 0.34 0.33 -0.68 0.42 1.00 1467 1537
## l1[2] -0.20 -0.20 0.24 0.25 -0.58 0.20 1.00 1078 1096
## l1[3] 0.24 0.24 0.19 0.20 -0.07 0.57 1.00 1194 1230
## l1[4] 1.63 1.63 0.25 0.26 1.20 2.04 1.00 1452 1564
## l1[5] -0.26 -0.26 0.25 0.26 -0.65 0.15 1.00 1076 1150
## l1[6] 0.95 0.95 0.23 0.24 0.56 1.33 1.00 1536 1334
## l1[7] 0.59 0.59 0.19 0.20 0.26 0.90 1.00 1397 1377
## l1[8] 1.01 1.01 0.24 0.25 0.61 1.40 1.00 1543 1359
## l1[9] 0.86 0.86 0.22 0.23 0.49 1.22 1.00 1517 1475
## l1[10] 0.07 0.08 0.30 0.30 -0.44 0.56 1.00 1509 1533
## l1[11] 0.13 0.14 0.29 0.29 -0.37 0.60 1.00 1525 1585
## l1[12] 0.00 0.00 0.22 0.22 -0.34 0.35 1.00 1101 1248
## l1[13] -0.67 -0.67 0.33 0.33 -1.19 -0.12 1.00 1082 1367
## l1[14] 0.04 0.05 0.31 0.30 -0.48 0.53 1.00 1499 1533
## l1[15] 0.20 0.20 0.20 0.20 -0.13 0.52 1.00 1166 1264
## l1[16] -0.17 -0.17 0.24 0.24 -0.54 0.22 1.00 1079 1122
## l1[17] -0.21 -0.20 0.36 0.35 -0.81 0.35 1.00 1452 1479
## l1[18] 0.50 0.51 0.24 0.24 0.10 0.88 1.00 1662 1639
## l1[19] 0.30 0.30 0.19 0.20 -0.01 0.62 1.00 1225 1230
## l1[20] 1.60 1.61 0.25 0.25 1.18 2.01 1.00 1464 1566
## l1[21] 0.01 0.01 0.21 0.22 -0.34 0.36 1.00 1102 1248
## l1[22] -0.93 -0.93 0.38 0.39 -1.54 -0.29 1.00 1099 1378
## l1[23] 0.20 0.20 0.20 0.20 -0.12 0.53 1.00 1166 1280
## l1[24] -0.80 -0.80 0.35 0.36 -1.36 -0.21 1.00 1091 1363
## l1[25] 0.82 0.82 0.21 0.22 0.46 1.17 1.00 1506 1545
## l1[26] 0.52 0.52 0.19 0.20 0.19 0.83 1.00 1358 1377
## l1[27] 0.57 0.58 0.23 0.23 0.19 0.94 1.00 1698 1654
## l1[28] 0.61 0.62 0.22 0.23 0.23 0.97 1.00 1716 1602
## l1[29] 1.38 1.39 0.22 0.23 1.01 1.75 1.00 1566 1383
## l1[30] -0.09 -0.09 0.23 0.23 -0.44 0.28 1.00 1086 1217
y1=mcmc$draws('y1')
smy=summary(y1)
table(tb$y,smy$median)##
## 0 1 2 3 4 5
## 0 2 6 0 0 0 0
## 1 1 5 1 0 0 0
## 2 0 5 1 1 1 0
## 3 0 0 4 0 0 0
## 4 0 0 0 0 0 1
## 6 0 0 1 0 0 1cat('\n')table(tb$y,smy$median,tb$c)## , , = a
##
##
## 0 1 2 3 4 5
## 0 2 4 0 0 0 0
## 1 1 4 1 0 0 0
## 2 0 2 1 1 0 0
## 3 0 0 2 0 0 0
## 4 0 0 0 0 0 0
## 6 0 0 1 0 0 0
##
## , , = b
##
##
## 0 1 2 3 4 5
## 0 0 2 0 0 0 0
## 1 0 1 0 0 0 0
## 2 0 3 0 0 1 0
## 3 0 0 2 0 0 0
## 4 0 0 0 0 0 1
## 6 0 0 0 0 0 1par(mfrow=c(1,2))
hist(tb$y-smy$median,xlab='obs.-prd.',main='residual')
density(tb$y-smy$median) |> plot(xlab='obs.-prd.',main='residual')
map=c()
for(i in 1:n){
a=table(y1[,,i])
map[i]=as.integer(names(a[a==max(a)]))
}
table(tb$y,map)## map
## 0 1 2 3 4
## 0 5 3 0 0 0
## 1 4 3 0 0 0
## 2 1 4 2 1 0
## 3 0 3 1 0 0
## 4 0 0 0 0 1
## 6 0 0 1 0 1cat('\n')table(tb$y,map,tb$c)## , , = a
##
## map
## 0 1 2 3 4
## 0 4 2 0 0 0
## 1 4 2 0 0 0
## 2 0 2 2 0 0
## 3 0 1 1 0 0
## 4 0 0 0 0 0
## 6 0 0 1 0 0
##
## , , = b
##
## map
## 0 1 2 3 4
## 0 1 1 0 0 0
## 1 0 1 0 0 0
## 2 1 2 0 1 0
## 3 0 2 0 0 0
## 4 0 0 0 0 1
## 6 0 0 0 0 1par(mfrow=c(1,2))
hist(tb$y-map,xlab='obs.-map',main='residual')
density(tb$y-map) |> plot(xlab='obs.-map',main='residual')
logistic regression
# of samples is N,
objective variable 0/1 binary
probability of incident pi[0,1]
logit pi=log(pi/ 1-pi)=sum(bj*xji),j=[0,K],i=[1,N] (-Infinity, Infinity)
yi~Ber(pi), 0/1 binary
odds(x)=p(x)/ 1-p(x), probablity of incident / probablity of no incident
odds ratio(x0,x1)=odds(x1)/odd(x0)
when xj -> xj +1, odds ratio -> odds ratio *exp bj ex6-2.stan
data{
int N;
int K;
array[N] int y;
matrix[N,K] X;
}
parameters{
vector[K] b;
}
model{
vector[N] z=X*b;
y~bernoulli_logit(z);
}
generated quantities{
array[N] int y1;
vector[N] z1=X*b;
for(i in 1:N){
y1[i]=bernoulli_rng(inv_logit(z1[i]));
}
}n=30
x=runif(n,-1,1)
c=sample(c('a','b'),n,replace=T)
z=x+(c=='b')
y=rbinom(n,1,1/(1+exp(-z)))
tb=tibble(x=x,c=c,y=y)
f=formula(y~x+c)
qplot(data=tb,x,y,col=c)
glm(f,tb,family='binomial') # it can caluculte when all trials are once##
## Call: glm(formula = f, family = "binomial", data = tb)
##
## Coefficients:
## (Intercept) x cb
## 0.018 4.034 2.779
##
## Degrees of Freedom: 29 Total (i.e. Null); 27 Residual
## Null Deviance: 32.6
## Residual Deviance: 20.7 AIC: 26.7X=model.matrix(f,tb)
data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X)
mdl=cmdstan_model('./ex6-2.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -28.9854
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 11 -10.3708 0.000328813 2.39543e-05 1 1 15
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -10.37
## b[1] 0.02
## b[2] 4.03
## b[3] 2.78
## y1[1] 0.00
## y1[2] 1.00
## y1[3] 1.00
## y1[4] 0.00
## y1[5] 1.00
## y1[6] 1.00
## y1[7] 0.00
## y1[8] 1.00
## y1[9] 1.00
## y1[10] 0.00
## y1[11] 1.00
## y1[12] 1.00
## y1[13] 0.00
## y1[14] 1.00
## y1[15] 1.00
## y1[16] 0.00
## y1[17] 1.00
## y1[18] 1.00
## y1[19] 1.00
## y1[20] 1.00
## y1[21] 1.00
## y1[22] 1.00
## y1[23] 1.00
## y1[24] 1.00
## y1[25] 1.00
## y1[26] 0.00
## y1[27] 1.00
## y1[28] 0.00
## y1[29] 1.00
## y1[30] 1.00
## z1[1] 1.68
## z1[2] 1.09
## z1[3] 2.46
## z1[4] -2.51
## z1[5] 0.61
## z1[6] 4.20
## z1[7] 0.61
## z1[8] 3.71
## z1[9] 2.52
## z1[10] 1.36
## z1[11] -0.69
## z1[12] 0.18
## z1[13] -1.02
## z1[14] 1.75
## z1[15] 3.27
## z1[16] 1.72
## z1[17] 6.08
## z1[18] 5.77
## z1[19] 1.09
## z1[20] 2.85
## z1[21] 4.88
## z1[22] 3.86
## z1[23] 6.13
## z1[24] 6.31
## z1[25] 4.12
## z1[26] 0.26
## z1[27] 0.80
## z1[28] -2.05
## z1[29] 2.03
## z1[30] 1.01mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc,exc='z1',ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -12.03 -11.67 1.35 1.07 -14.84 -10.58 1.00 666 807
## b[1] -0.13 -0.08 0.96 0.94 -1.81 1.35 1.01 710 759
## b[2] 5.28 5.00 2.19 2.07 2.19 9.30 1.00 414 544
## b[3] 3.72 3.50 1.95 1.83 0.95 7.34 1.01 397 544
## y1[1] 0.86 1.00 0.35 0.00 0.00 1.00 1.00 1908 NA
## y1[2] 0.75 1.00 0.43 0.00 0.00 1.00 1.00 1934 NA
## y1[3] 0.92 1.00 0.27 0.00 0.00 1.00 1.00 1980 NA
## y1[4] 0.09 0.00 0.28 0.00 0.00 1.00 1.00 1131 NA
## y1[5] 0.65 1.00 0.48 0.00 0.00 1.00 1.00 1829 NA
## y1[6] 0.98 1.00 0.15 0.00 1.00 1.00 1.00 1849 NA
## y1[7] 0.66 1.00 0.47 0.00 0.00 1.00 1.00 1761 NA
## y1[8] 0.98 1.00 0.15 0.00 1.00 1.00 1.00 1772 NA
## y1[9] 0.93 1.00 0.26 0.00 0.00 1.00 1.00 1723 NA
## y1[10] 0.81 1.00 0.39 0.00 0.00 1.00 1.00 1881 NA
## y1[11] 0.31 0.00 0.46 0.00 0.00 1.00 1.00 1924 NA
## y1[12] 0.55 1.00 0.50 0.00 0.00 1.00 1.00 2037 NA
## y1[13] 0.24 0.00 0.42 0.00 0.00 1.00 1.00 1332 NA
## y1[14] 0.86 1.00 0.35 0.00 0.00 1.00 1.00 1783 NA
## y1[15] 0.96 1.00 0.20 0.00 1.00 1.00 1.00 1938 NA
## y1[16] 0.86 1.00 0.34 0.00 0.00 1.00 1.00 1931 NA
## y1[17] 1.00 1.00 0.07 0.00 1.00 1.00 1.00 1689 NA
## y1[18] 1.00 1.00 0.06 0.00 1.00 1.00 1.00 2030 NA
## y1[19] 0.77 1.00 0.42 0.00 0.00 1.00 1.00 1961 NA
## y1[20] 0.95 1.00 0.22 0.00 1.00 1.00 1.00 1835 NA
## y1[21] 0.99 1.00 0.10 0.00 1.00 1.00 1.00 1717 NA
## y1[22] 0.98 1.00 0.15 0.00 1.00 1.00 1.00 1712 NA
## y1[23] 0.99 1.00 0.07 0.00 1.00 1.00 1.00 1192 NA
## y1[24] 1.00 1.00 0.07 0.00 1.00 1.00 1.00 1736 NA
## y1[25] 0.98 1.00 0.13 0.00 1.00 1.00 1.00 1932 NA
## y1[26] 0.56 1.00 0.50 0.00 0.00 1.00 1.00 2165 NA
## y1[27] 0.70 1.00 0.46 0.00 0.00 1.00 1.00 2141 NA
## y1[28] 0.13 0.00 0.33 0.00 0.00 1.00 1.00 1548 NA
## y1[29] 0.88 1.00 0.32 0.00 0.00 1.00 1.00 2003 NA
## y1[30] 0.76 1.00 0.43 0.00 0.00 1.00 1.00 1941 NA
## z1[1] 2.04 1.95 1.01 0.95 0.57 3.78 1.00 1098 1137
## z1[2] 1.27 1.21 0.89 0.87 -0.07 2.78 1.00 1201 1203
## z1[3] 3.07 2.97 1.27 1.23 1.24 5.29 1.00 771 872
## z1[4] -3.45 -3.20 1.98 1.87 -7.05 -0.65 1.01 411 544
## z1[5] 0.64 0.62 0.88 0.86 -0.75 2.09 1.00 1018 1117
## z1[6] 5.43 5.15 2.07 1.90 2.58 9.28 1.00 392 559
## z1[7] 0.73 0.70 0.84 0.83 -0.65 2.15 1.00 1878 1404
## z1[8] 4.70 4.53 1.83 1.78 2.10 7.97 1.00 578 761
## z1[9] 3.22 3.07 1.29 1.20 1.41 5.55 1.00 455 651
## z1[10] 1.70 1.64 0.91 0.89 0.29 3.29 1.00 806 979
## z1[11] -0.98 -0.94 1.14 1.08 -2.93 0.77 1.00 1055 1055
## z1[12] 0.16 0.14 0.89 0.86 -1.28 1.61 1.00 1891 1443
## z1[13] -1.49 -1.33 1.31 1.20 -3.83 0.43 1.01 482 607
## z1[14] 2.14 2.05 1.03 0.97 0.64 3.91 1.00 1064 1137
## z1[15] 4.21 3.99 1.62 1.49 1.96 7.20 1.00 410 569
## z1[16] 2.18 2.08 1.00 0.97 0.70 3.95 1.00 611 764
## z1[17] 7.89 7.48 3.03 2.78 3.76 13.46 1.00 381 541
## z1[18] 7.48 7.12 2.87 2.64 3.55 12.74 1.00 382 541
## z1[19] 1.27 1.21 0.89 0.87 -0.07 2.78 1.00 1201 1203
## z1[20] 3.57 3.44 1.43 1.39 1.54 6.05 1.00 683 711
## z1[21] 6.32 6.00 2.41 2.20 3.00 10.84 1.00 385 556
## z1[22] 4.89 4.72 1.90 1.85 2.20 8.30 1.00 567 706
## z1[23] 7.95 7.52 3.05 2.80 3.77 13.56 1.00 381 541
## z1[24] 8.19 7.74 3.15 2.87 3.88 13.99 1.00 381 541
## z1[25] 5.32 5.05 2.03 1.86 2.52 9.08 1.00 393 559
## z1[26] 0.27 0.25 0.88 0.84 -1.14 1.72 1.00 1953 1519
## z1[27] 0.97 0.94 0.84 0.82 -0.41 2.39 1.00 1561 1354
## z1[28] -2.84 -2.63 1.76 1.66 -5.96 -0.34 1.01 423 533
## z1[29] 2.51 2.42 1.12 1.07 0.90 4.40 1.00 920 984
## z1[30] 1.25 1.21 0.86 0.83 -0.11 2.71 1.00 1185 1251
y1=mcmc$draws('y1')
smy=summary(y1)
table(tb$y,smy$median)##
## 0 1
## 0 3 4
## 1 1 22cat('\n')table(tb$y,smy$median,tb$c)## , , = a
##
##
## 0 1
## 0 3 1
## 1 0 9
##
## , , = b
##
##
## 0 1
## 0 0 3
## 1 1 13par(mfrow=c(1,2))
hist(tb$y-smy$median,xlab='obs.-prd.',main='residual')
density(tb$y-smy$median) |> plot(xlab='obs.-prd.',main='residual')
map=c()
for(i in 1:n){
a=table(y1[,,i])
map[i]=as.integer(names(a[a==max(a)]))
}
table(tb$y,map)## map
## 0 1
## 0 3 4
## 1 1 22cat('\n')table(tb$y,map,tb$c)## , , = a
##
## map
## 0 1
## 0 3 1
## 1 0 9
##
## , , = b
##
## map
## 0 1
## 0 0 3
## 1 1 13par(mfrow=c(1,2))
hist(tb$y-map,xlab='obs.-map',main='residual')
density(tb$y-map) |> plot(xlab='obs.-map',main='residual')
multi logistic regression
# of samples is N,
# of trials of a sample i is mi,
objective variable [0,n], integer
probability of incident pi[0,1]
logit pi=log(pi/ 1-pi)=sum(bj*xji),j=[0,K],i=[1,N] (-Infinity, Infinity)
yi~B(mi, pi), # of acutual incident
odds(x)=p(x)/ 1-p(x), probablity of incident / probablity of no incident
odds ratio(x0,x1)=odds(x1)/odd(x0)
when xj -> xj +1, odds ratio -> odds ratio *exp bj ex6-3.stan
data{
int N;
int K;
array[N] int m;
array[N] int y;
matrix[N,K] X;
}
parameters{
vector[K] b;
}
model{
vector[N] z=X*b;
y~binomial_logit(m,z);
}
generated quantities{
array[N] int y1;
vector[N] z1=X*b;
for(i in 1:N){
y1[i]=binomial_rng(m[i],inv_logit(z1[i]));
}
}n=30
m=floor(runif(n,1,10)) # trials are varying (1,10)
x=runif(n,-1,1)
c=sample(c('a','b'),n,replace=T)
z=x+(c=='b')
y=rbinom(n,m,1/(1+exp(-z)))
tb=tibble(x=x,c=c,y=y)
f=formula(y~x+c)
qplot(data=tb,x,y,col=c)
X=model.matrix(f,tb)
data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X,m=m)
mdl=cmdstan_model('./ex6-3.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -180.062
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 9 -114.455 0.00117686 0.00149793 1 1 11
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -114.46
## b[1] 0.12
## b[2] 0.43
## b[3] 0.38
## y1[1] 4.00
## y1[2] 3.00
## y1[3] 2.00
## y1[4] 4.00
## y1[5] 0.00
## y1[6] 2.00
## y1[7] 6.00
## y1[8] 5.00
## y1[9] 4.00
## y1[10] 4.00
## y1[11] 1.00
## y1[12] 1.00
## y1[13] 4.00
## y1[14] 2.00
## y1[15] 4.00
## y1[16] 4.00
## y1[17] 4.00
## y1[18] 6.00
## y1[19] 4.00
## y1[20] 7.00
## y1[21] 4.00
## y1[22] 5.00
## y1[23] 2.00
## y1[24] 2.00
## y1[25] 8.00
## y1[26] 3.00
## y1[27] 2.00
## y1[28] 4.00
## y1[29] 0.00
## y1[30] 1.00
## z1[1] 0.43
## z1[2] 0.09
## z1[3] 0.15
## z1[4] 0.36
## z1[5] -0.06
## z1[6] -0.27
## z1[7] -0.23
## z1[8] 0.77
## z1[9] -0.11
## z1[10] 0.56
## z1[11] 0.04
## z1[12] 0.59
## z1[13] 0.87
## z1[14] 0.25
## z1[15] 0.08
## z1[16] 0.34
## z1[17] 0.91
## z1[18] 0.23
## z1[19] 0.45
## z1[20] 0.33
## z1[21] 0.42
## z1[22] 0.40
## z1[23] 0.47
## z1[24] 0.71
## z1[25] 0.87
## z1[26] 0.53
## z1[27] 0.42
## z1[28] 0.61
## z1[29] -0.04
## z1[30] 0.22mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc,exc='z1',ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -116.01 -115.68 1.33 1.00 -118.54 -114.62 1.00 826 1232
## b[1] 0.12 0.13 0.22 0.21 -0.24 0.48 1.00 1264 1200
## b[2] 0.44 0.45 0.28 0.29 -0.02 0.89 1.00 1353 1256
## b[3] 0.38 0.38 0.32 0.32 -0.12 0.91 1.00 972 1007
## y1[1] 3.61 4.00 1.24 1.48 2.00 6.00 1.00 2026 NA
## y1[2] 3.06 3.00 1.27 1.48 1.00 5.00 1.00 1929 1998
## y1[3] 1.06 1.00 0.71 1.48 0.00 2.00 1.00 1956 NA
## y1[4] 2.34 2.00 0.99 1.48 1.00 4.00 1.00 1853 NA
## y1[5] 2.94 3.00 1.26 1.48 1.00 5.00 1.00 1849 1973
## y1[6] 2.56 3.00 1.30 1.48 0.95 5.00 1.00 2094 2053
## y1[7] 3.49 3.00 1.52 1.48 1.00 6.00 1.00 1806 1836
## y1[8] 5.45 6.00 1.41 1.48 3.00 8.00 1.00 1836 NA
## y1[9] 2.83 3.00 1.27 1.48 1.00 5.00 1.00 1857 1890
## y1[10] 3.18 3.00 1.12 1.48 1.00 5.00 1.00 1913 NA
## y1[11] 1.56 2.00 0.88 1.48 0.00 3.00 1.00 1985 NA
## y1[12] 1.28 1.00 0.68 1.48 0.00 2.00 1.00 1853 NA
## y1[13] 4.22 4.00 1.17 1.48 2.00 6.00 1.00 1961 NA
## y1[14] 4.53 5.00 1.44 1.48 2.00 7.00 1.00 1785 1893
## y1[15] 2.08 2.00 1.05 1.48 0.00 4.00 1.00 1863 NA
## y1[16] 4.73 5.00 1.46 1.48 2.00 7.00 1.00 1926 1953
## y1[17] 3.58 4.00 1.06 1.48 2.00 5.00 1.00 1774 NA
## y1[18] 5.04 5.00 1.66 1.48 2.00 8.00 1.00 1724 1743
## y1[19] 5.52 6.00 1.54 1.48 3.00 8.00 1.00 1990 1796
## y1[20] 5.22 5.00 1.57 1.48 3.00 8.00 1.00 1905 1683
## y1[21] 3.04 3.00 1.11 1.48 1.00 5.00 1.00 1758 NA
## y1[22] 4.19 4.00 1.38 1.48 2.00 6.00 1.00 2076 1956
## y1[23] 2.49 3.00 0.98 1.48 1.00 4.00 1.00 1924 NA
## y1[24] 2.69 3.00 0.96 1.48 1.00 4.00 1.00 1713 NA
## y1[25] 6.35 6.00 1.42 1.48 4.00 9.00 1.00 1817 NA
## y1[26] 3.78 4.00 1.25 1.48 2.00 6.00 1.00 1833 NA
## y1[27] 1.21 1.00 0.69 1.48 0.00 2.00 1.00 1608 NA
## y1[28] 5.90 6.00 1.48 1.48 3.00 8.00 1.00 1784 1784
## y1[29] 1.48 1.00 0.87 1.48 0.00 3.00 1.00 2088 NA
## y1[30] 1.66 2.00 0.89 1.48 0.00 3.00 1.00 1947 NA
## z1[1] 0.44 0.43 0.27 0.27 0.00 0.89 1.00 1369 1538
## z1[2] 0.09 0.10 0.22 0.21 -0.27 0.45 1.00 1259 1333
## z1[3] 0.16 0.15 0.36 0.36 -0.42 0.74 1.00 1427 1259
## z1[4] 0.37 0.37 0.27 0.28 -0.07 0.81 1.00 1527 1618
## z1[5] -0.06 -0.05 0.26 0.25 -0.50 0.37 1.00 1252 1127
## z1[6] -0.27 -0.27 0.35 0.34 -0.85 0.31 1.00 1279 1233
## z1[7] -0.23 -0.23 0.33 0.32 -0.77 0.32 1.00 1271 1165
## z1[8] 0.79 0.78 0.27 0.26 0.37 1.24 1.00 1555 1744
## z1[9] -0.11 -0.11 0.28 0.27 -0.58 0.34 1.00 1255 1138
## z1[10] 0.57 0.57 0.23 0.23 0.20 0.96 1.00 1578 1616
## z1[11] 0.04 0.05 0.23 0.22 -0.34 0.42 1.00 1252 1129
## z1[12] 0.60 0.60 0.23 0.23 0.23 0.98 1.00 1572 1692
## z1[13] 0.89 0.88 0.30 0.31 0.40 1.40 1.00 1496 1672
## z1[14] 0.26 0.26 0.22 0.21 -0.10 0.62 1.00 1295 1264
## z1[15] 0.08 0.08 0.40 0.40 -0.56 0.73 1.00 1410 1276
## z1[16] 0.35 0.35 0.24 0.24 -0.04 0.75 1.00 1333 1426
## z1[17] 0.93 0.92 0.32 0.33 0.42 1.48 1.00 1492 1586
## z1[18] 0.24 0.24 0.32 0.33 -0.29 0.76 1.00 1453 1273
## z1[19] 0.46 0.45 0.28 0.28 0.00 0.93 1.00 1375 1602
## z1[20] 0.34 0.33 0.24 0.23 -0.05 0.73 1.00 1329 1403
## z1[21] 0.43 0.42 0.27 0.26 -0.01 0.87 1.00 1366 1591
## z1[22] 0.41 0.41 0.26 0.26 -0.01 0.83 1.00 1556 1526
## z1[23] 0.48 0.47 0.24 0.24 0.08 0.87 1.00 1604 1457
## z1[24] 0.73 0.72 0.25 0.25 0.33 1.15 1.00 1570 1781
## z1[25] 0.89 0.88 0.31 0.31 0.40 1.41 1.00 1494 1672
## z1[26] 0.54 0.53 0.32 0.31 0.04 1.07 1.00 1392 1598
## z1[27] 0.43 0.43 0.25 0.25 0.02 0.84 1.00 1569 1520
## z1[28] 0.63 0.62 0.23 0.23 0.25 1.02 1.00 1567 1739
## z1[29] -0.04 -0.03 0.25 0.24 -0.46 0.38 1.00 1253 1062
## z1[30] 0.22 0.22 0.22 0.21 -0.13 0.57 1.00 1283 1197
y1=mcmc$draws('y1')
smy=summary(y1)
table(tb$y,smy$median)##
## 1 2 3 4 5 6
## 1 2 2 0 0 0 0
## 2 2 1 2 1 0 0
## 3 0 0 5 1 1 0
## 4 0 1 2 2 0 0
## 5 0 0 0 1 2 2
## 6 0 0 0 0 1 1
## 7 0 0 0 0 0 1cat('\n')table(tb$y,smy$median,tb$c)## , , = a
##
##
## 1 2 3 4 5 6
## 1 0 2 0 0 0 0
## 2 1 0 1 0 0 0
## 3 0 0 3 0 1 0
## 4 0 0 2 2 0 0
## 5 0 0 0 0 2 0
## 6 0 0 0 0 0 1
## 7 0 0 0 0 0 0
##
## , , = b
##
##
## 1 2 3 4 5 6
## 1 2 0 0 0 0 0
## 2 1 1 1 1 0 0
## 3 0 0 2 1 0 0
## 4 0 1 0 0 0 0
## 5 0 0 0 1 0 2
## 6 0 0 0 0 1 0
## 7 0 0 0 0 0 1par(mfrow=c(1,2))
hist(tb$y-smy$median,xlab='obs.-prd.',main='residual')
density(tb$y-smy$median) |> plot(xlab='obs.-prd.',main='residual')
map=c()
for(i in 1:n){
a=table(y1[,,i])
map[i]=as.integer(names(a[a==max(a)]))
}
table(tb$y,map)## map
## 1 2 3 4 5 6 7
## 1 2 2 0 0 0 0 0
## 2 2 1 2 1 0 0 0
## 3 0 0 5 2 0 0 0
## 4 0 1 2 2 0 0 0
## 5 0 0 0 0 3 2 0
## 6 0 0 0 0 1 1 0
## 7 0 0 0 0 0 0 1cat('\n')table(tb$y,map,tb$c)## , , = a
##
## map
## 1 2 3 4 5 6 7
## 1 0 2 0 0 0 0 0
## 2 1 0 1 0 0 0 0
## 3 0 0 3 1 0 0 0
## 4 0 0 2 2 0 0 0
## 5 0 0 0 0 2 0 0
## 6 0 0 0 0 0 1 0
## 7 0 0 0 0 0 0 0
##
## , , = b
##
## map
## 1 2 3 4 5 6 7
## 1 2 0 0 0 0 0 0
## 2 1 1 1 1 0 0 0
## 3 0 0 2 1 0 0 0
## 4 0 1 0 0 0 0 0
## 5 0 0 0 0 1 2 0
## 6 0 0 0 0 1 0 0
## 7 0 0 0 0 0 0 1par(mfrow=c(1,2))
hist(tb$y-map,xlab='obs.-map',main='residual')
density(tb$y-map) |> plot(xlab='obs.-map',main='residual')
ex7 interaction of variables
n=50
mdl=cmdstan_model('./ex5.stan')
tb=tibble(x=runif(n,-3,3),
ca=sample(c('a','b'),n,replace=T),cb=sample(c('a','b'),n,replace=T),
y=rnorm(n,x+(ca=='b')+(cb=='b')-(ca=='b')*(cb=='b'),1))
grid.arrange(qplot(data=tb,x,y,col=ca),
qplot(data=tb,x,y,col=cb),ncol=2)
fn=function(f){
X=model.matrix(f,tb)
data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X)
mle=mdl$optimize(data=data)
mle
mcmc=goMCMC(mdl,data)
seeMCMC(mcmc,exc='m',ch=F)
y1=mcmc$draws('y1')
smy=summary(y1)
tb$y1=smy$median
tb$l5=smy$q5
tb$u5=smy$q95
grid.arrange(
qplot(data=tb,y,y1,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))
par(mfrow=c(1,2))
hist(tb$y-tb$y1,xlab='obs.-prd.',main='residual')
density(tb$y-tb$y1) |> plot(xlab='obs.-prd.',main='residual')
tb=arrange(tb,y1)
grid.arrange(
qplot(data=tb,1:n,y,col=ca)+
geom_point(aes(x=1:n,y=y1),tb,col='black')+
geom_line(aes(x=1:n,y=l5),tb,col='gray')+
geom_line(aes(x=1:n,y=u5),tb,col='gray'),
qplot(data=tb,1:n,y,col=cb)+
geom_point(aes(x=1:n,y=y1),tb,col='black')+
geom_line(aes(x=1:n,y=l5),tb,col='gray')+
geom_line(aes(x=1:n,y=u5),tb,col='gray'),
ncol=2
)
}
f0=formula(y~x+ca)
f1=formula(y~x+ca+cb)
f2=formula(y~x+ca*cb)
fn(f0)## Initial log joint probability = -613.308
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 14 -30.2492 0.00026474 0.000991645 0.9285 0.9285 15
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.
## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
##
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -32.27 -31.94 1.50 1.24 -35.31 -30.48 1.00 918 1472
## b[1] 0.34 0.34 0.23 0.22 -0.06 0.70 1.01 755 948
## b[2] 0.90 0.90 0.11 0.11 0.72 1.08 1.00 2517 1595
## b[3] 0.68 0.68 0.34 0.34 0.12 1.26 1.01 608 707
## s 1.18 1.17 0.13 0.12 0.99 1.41 1.00 2424 1683
## y1[1] -0.76 -0.73 1.22 1.22 -2.84 1.18 1.00 1904 1933
## y1[2] 0.15 0.11 1.24 1.27 -1.81 2.20 1.00 1743 1599
## y1[3] 0.43 0.42 1.21 1.19 -1.57 2.46 1.00 1789 2011
## y1[4] 0.58 0.58 1.18 1.18 -1.35 2.47 1.00 1803 2054
## y1[5] -0.29 -0.30 1.22 1.21 -2.24 1.72 1.00 1771 1704
## y1[6] -0.16 -0.16 1.22 1.17 -2.22 1.85 1.00 2031 2016
## y1[7] 0.18 0.20 1.26 1.22 -1.89 2.22 1.00 1783 1856
## y1[8] 3.38 3.39 1.24 1.25 1.34 5.53 1.00 2173 2055
## y1[9] -1.88 -1.90 1.24 1.18 -3.89 0.14 1.00 2059 1970
## y1[10] 0.35 0.36 1.21 1.20 -1.65 2.33 1.00 1914 1893
## y1[11] 0.47 0.47 1.22 1.20 -1.50 2.49 1.00 1869 1792
## y1[12] 1.16 1.17 1.18 1.13 -0.82 3.15 1.00 2006 1753
## y1[13] 2.08 2.09 1.22 1.22 0.08 4.07 1.00 2063 2129
## y1[14] 1.77 1.79 1.23 1.19 -0.25 3.74 1.00 1986 1966
## y1[15] 1.91 1.89 1.24 1.17 -0.12 3.97 1.00 1854 1865
## y1[16] 3.06 3.03 1.23 1.20 1.09 5.08 1.00 1953 1928
## y1[17] 3.00 2.98 1.24 1.25 1.03 5.04 1.00 1991 1852
## y1[18] 1.58 1.58 1.24 1.20 -0.43 3.66 1.00 2020 1657
## y1[19] -0.33 -0.33 1.24 1.21 -2.36 1.72 1.00 1944 2057
## y1[20] 3.53 3.51 1.24 1.20 1.39 5.57 1.00 2117 1955
## y1[21] -1.59 -1.64 1.19 1.18 -3.52 0.37 1.00 1820 1973
## y1[22] 0.23 0.25 1.17 1.13 -1.67 2.15 1.00 1941 1712
## y1[23] 1.95 1.96 1.21 1.19 -0.06 3.90 1.00 1577 1527
## y1[24] 2.61 2.62 1.24 1.27 0.56 4.66 1.00 1558 1694
## y1[25] 1.09 1.10 1.24 1.17 -1.00 3.08 1.00 2028 1914
## y1[26] -1.38 -1.39 1.20 1.25 -3.29 0.58 1.00 2065 1967
## y1[27] -0.10 -0.03 1.23 1.21 -2.16 1.88 1.00 1994 1831
## y1[28] 2.13 2.12 1.19 1.18 0.19 4.04 1.00 2023 2013
## y1[29] 2.31 2.32 1.21 1.21 0.39 4.38 1.00 1912 1779
## y1[30] 1.70 1.70 1.19 1.20 -0.33 3.63 1.00 1960 1971
## y1[31] -1.11 -1.11 1.23 1.16 -3.16 0.93 1.00 1857 1811
## y1[32] -1.50 -1.51 1.23 1.20 -3.53 0.53 1.00 2004 1933
## y1[33] 2.55 2.52 1.28 1.24 0.44 4.67 1.00 1746 1841
## y1[34] 0.21 0.23 1.21 1.24 -1.82 2.25 1.00 1915 1630
## y1[35] 2.67 2.66 1.19 1.15 0.66 4.63 1.00 1771 1822
## y1[36] -1.35 -1.35 1.23 1.18 -3.44 0.65 1.00 1847 1781
## y1[37] -0.28 -0.31 1.23 1.26 -2.32 1.66 1.00 1827 1672
## y1[38] 0.35 0.35 1.19 1.16 -1.60 2.30 1.00 1828 1855
## y1[39] -2.18 -2.19 1.25 1.18 -4.21 -0.16 1.00 1782 1623
## y1[40] -0.29 -0.27 1.19 1.21 -2.24 1.62 1.00 1863 1858
## y1[41] 0.89 0.94 1.22 1.18 -1.09 2.91 1.00 1968 1892
## y1[42] -0.95 -0.94 1.26 1.24 -3.06 1.11 1.00 1862 1933
## y1[43] 3.10 3.10 1.23 1.21 1.10 5.13 1.00 2015 1894
## y1[44] 1.42 1.44 1.21 1.17 -0.60 3.44 1.00 1905 1928
## y1[45] 1.91 1.93 1.18 1.17 -0.02 3.90 1.00 2188 1664
## y1[46] -0.73 -0.76 1.24 1.21 -2.75 1.32 1.00 2122 1839
## y1[47] 0.00 -0.03 1.22 1.16 -1.97 2.12 1.00 1727 1868
## y1[48] 2.04 2.01 1.21 1.15 0.08 4.06 1.00 1876 1954
## y1[49] 1.81 1.82 1.19 1.14 -0.22 3.67 1.00 1832 1884
## y1[50] 2.92 2.92 1.20 1.22 0.94 4.85 1.00 2095 1973
## m1[1] -0.71 -0.71 0.26 0.26 -1.14 -0.28 1.00 996 1185
## m1[2] 0.11 0.12 0.23 0.22 -0.28 0.48 1.01 778 950
## m1[3] 0.41 0.42 0.23 0.22 0.02 0.78 1.01 751 988
## m1[4] 0.56 0.57 0.23 0.23 0.16 0.93 1.01 750 901
## m1[5] -0.34 -0.34 0.33 0.33 -0.86 0.19 1.00 1189 1411
## m1[6] -0.14 -0.13 0.23 0.23 -0.53 0.24 1.01 824 977
## m1[7] 0.23 0.23 0.29 0.29 -0.22 0.70 1.00 1081 1309
## m1[8] 3.33 3.33 0.34 0.33 2.77 3.89 1.00 1762 1532
## m1[9] -1.88 -1.88 0.35 0.36 -2.45 -1.33 1.00 1525 1540
## m1[10] 0.35 0.35 0.28 0.28 -0.09 0.81 1.00 1063 1327
## m1[11] 0.53 0.53 0.23 0.23 0.13 0.89 1.01 750 913
## m1[12] 1.17 1.17 0.25 0.25 0.76 1.57 1.00 1022 1448
## m1[13] 2.09 2.09 0.26 0.25 1.68 2.52 1.00 1221 1593
## m1[14] 1.76 1.76 0.25 0.24 1.35 2.17 1.00 1115 1426
## m1[15] 1.88 1.88 0.30 0.30 1.39 2.37 1.00 936 1162
## m1[16] 3.06 3.05 0.32 0.31 2.53 3.59 1.00 1633 1522
## m1[17] 2.97 2.96 0.31 0.30 2.46 3.49 1.00 1587 1558
## m1[18] 1.60 1.60 0.25 0.24 1.19 2.01 1.00 1078 1432
## m1[19] -0.38 -0.38 0.33 0.33 -0.90 0.17 1.00 1196 1390
## m1[20] 3.56 3.55 0.36 0.35 2.97 4.15 1.00 1860 1529
## m1[21] -1.53 -1.53 0.32 0.32 -2.04 -1.02 1.00 1359 1466
## m1[22] 0.23 0.24 0.23 0.22 -0.16 0.60 1.01 763 968
## m1[23] 1.93 1.94 0.31 0.31 1.44 2.43 1.00 949 1191
## m1[24] 2.58 2.58 0.37 0.37 1.99 3.18 1.00 1104 1328
## m1[25] 1.13 1.14 0.25 0.25 0.72 1.55 1.00 796 972
## m1[26] -1.39 -1.39 0.31 0.31 -1.89 -0.90 1.00 1292 1451
## m1[27] -0.07 -0.07 0.23 0.23 -0.46 0.30 1.01 811 987
## m1[28] 2.13 2.13 0.26 0.25 1.71 2.56 1.00 1234 1619
## m1[29] 2.32 2.31 0.27 0.26 1.88 2.76 1.00 1308 1549
## m1[30] 1.72 1.73 0.29 0.29 1.25 2.19 1.00 901 1162
## m1[31] -1.13 -1.13 0.29 0.29 -1.61 -0.66 1.00 1172 1350
## m1[32] -1.51 -1.51 0.44 0.44 -2.21 -0.80 1.00 1440 1496
## m1[33] 2.48 2.48 0.36 0.36 1.90 3.06 1.00 1082 1227
## m1[34] 0.23 0.23 0.29 0.29 -0.22 0.70 1.00 1081 1309
## m1[35] 2.67 2.67 0.29 0.28 2.20 3.15 1.00 1461 1545
## m1[36] -1.34 -1.34 0.42 0.42 -2.01 -0.66 1.00 1407 1454
## m1[37] -0.27 -0.26 0.24 0.24 -0.66 0.12 1.01 855 1030
## m1[38] 0.35 0.35 0.23 0.22 -0.05 0.71 1.01 754 948
## m1[39] -2.16 -2.17 0.38 0.39 -2.78 -1.56 1.00 1651 1495
## m1[40] -0.30 -0.29 0.24 0.24 -0.70 0.09 1.01 864 1022
## m1[41] 0.94 0.94 0.24 0.24 0.53 1.33 1.00 771 924
## m1[42] -0.91 -0.92 0.38 0.38 -1.52 -0.30 1.00 1312 1442
## m1[43] 3.13 3.13 0.32 0.32 2.59 3.67 1.00 1664 1484
## m1[44] 1.44 1.44 0.25 0.24 1.03 1.85 1.00 1050 1422
## m1[45] 1.90 1.90 0.31 0.30 1.40 2.39 1.00 941 1149
## m1[46] -0.76 -0.75 0.26 0.26 -1.19 -0.32 1.00 1014 1215
## m1[47] 0.04 0.05 0.23 0.22 -0.35 0.41 1.01 789 977
## m1[48] 2.03 2.03 0.26 0.25 1.62 2.45 1.00 1198 1593
## m1[49] 1.82 1.82 0.25 0.24 1.41 2.23 1.00 1131 1427
## m1[50] 2.93 2.93 0.31 0.30 2.43 3.45 1.00 1570 1558
fn(f1)## Initial log joint probability = -85.83
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 16 -26.1273 0.000562122 0.00136671 1 1 27
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.
## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
##
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -28.81 -28.44 1.69 1.45 -32.08 -26.73 1.00 878 927
## b[1] -0.04 -0.03 0.26 0.25 -0.48 0.40 1.00 1038 1177
## b[2] 0.92 0.92 0.10 0.10 0.75 1.09 1.00 2389 1402
## b[3] 0.46 0.45 0.32 0.31 -0.06 1.00 1.00 844 1257
## b[4] 0.91 0.91 0.32 0.31 0.38 1.42 1.00 1099 1149
## s 1.10 1.08 0.12 0.12 0.92 1.31 1.00 2119 1463
## y1[1] -0.16 -0.16 1.14 1.11 -2.04 1.68 1.00 1757 1967
## y1[2] -0.27 -0.32 1.16 1.08 -2.12 1.66 1.00 1881 1475
## y1[3] 0.98 1.02 1.14 1.13 -0.87 2.84 1.00 2000 1810
## y1[4] 0.16 0.17 1.12 1.09 -1.69 1.99 1.00 1844 1895
## y1[5] -0.05 -0.04 1.18 1.16 -1.97 1.88 1.00 2007 1846
## y1[6] -0.55 -0.58 1.14 1.14 -2.38 1.30 1.00 1942 1932
## y1[7] 0.52 0.51 1.09 1.07 -1.20 2.34 1.00 2033 1972
## y1[8] 2.81 2.82 1.16 1.16 0.88 4.74 1.00 1914 1921
## y1[9] -1.38 -1.38 1.15 1.15 -3.29 0.46 1.00 2010 1925
## y1[10] -0.29 -0.30 1.18 1.17 -2.27 1.66 1.00 1972 1678
## y1[11] 0.15 0.16 1.12 1.09 -1.64 1.96 1.00 1886 1805
## y1[12] 1.50 1.52 1.11 1.16 -0.38 3.28 1.00 2099 2013
## y1[13] 2.46 2.48 1.15 1.11 0.57 4.35 1.00 1939 1715
## y1[14] 1.18 1.16 1.12 1.05 -0.68 2.99 1.00 1962 1973
## y1[15] 2.48 2.46 1.18 1.15 0.58 4.44 1.00 1831 1819
## y1[16] 3.43 3.43 1.15 1.14 1.59 5.38 1.00 1888 1874
## y1[17] 2.44 2.46 1.13 1.10 0.56 4.24 1.00 2041 1933
## y1[18] 0.97 1.02 1.16 1.14 -0.93 2.86 1.01 1876 1886
## y1[19] -0.13 -0.11 1.17 1.14 -2.05 1.76 1.00 1862 1722
## y1[20] 3.91 3.89 1.15 1.08 2.08 5.85 1.00 2010 1798
## y1[21] -1.92 -1.90 1.15 1.11 -3.86 -0.01 1.00 1803 1695
## y1[22] 0.84 0.84 1.12 1.10 -1.08 2.72 1.00 1870 1783
## y1[23] 1.59 1.55 1.15 1.12 -0.26 3.53 1.00 1860 1890
## y1[24] 3.14 3.14 1.19 1.20 1.17 5.01 1.00 1891 1863
## y1[25] 1.68 1.69 1.15 1.13 -0.22 3.57 1.00 1907 1833
## y1[26] -0.89 -0.89 1.14 1.15 -2.73 0.90 1.00 1899 1737
## y1[27] -0.45 -0.45 1.15 1.13 -2.30 1.45 1.00 1853 1882
## y1[28] 2.48 2.48 1.15 1.14 0.59 4.31 1.00 1926 1918
## y1[29] 2.65 2.63 1.14 1.10 0.81 4.57 1.00 1908 1768
## y1[30] 2.25 2.26 1.15 1.13 0.41 4.17 1.00 2060 1695
## y1[31] -1.59 -1.57 1.14 1.13 -3.51 0.29 1.00 1918 1961
## y1[32] -2.19 -2.21 1.18 1.15 -4.16 -0.28 1.00 1974 1784
## y1[33] 2.17 2.20 1.17 1.14 0.21 4.04 1.00 1922 1932
## y1[34] 0.52 0.52 1.12 1.14 -1.30 2.34 1.00 1818 1834
## y1[35] 2.12 2.09 1.16 1.16 0.29 4.01 1.00 1845 1818
## y1[36] -1.06 -1.07 1.19 1.16 -2.99 0.88 1.00 1873 1882
## y1[37] -0.70 -0.69 1.14 1.14 -2.55 1.16 1.00 1939 1850
## y1[38] -0.02 0.00 1.14 1.12 -1.94 1.76 1.00 1804 1956
## y1[39] -2.61 -2.61 1.18 1.13 -4.55 -0.69 1.00 1942 1813
## y1[40] 0.21 0.22 1.13 1.10 -1.69 2.01 1.00 1770 1921
## y1[41] 0.62 0.60 1.11 1.07 -1.22 2.43 1.00 1858 1869
## y1[42] -0.67 -0.69 1.17 1.16 -2.59 1.23 1.00 1884 1937
## y1[43] 2.59 2.57 1.14 1.15 0.75 4.52 1.00 1925 2013
## y1[44] 1.76 1.76 1.13 1.10 -0.04 3.66 1.00 1920 2014
## y1[45] 2.47 2.43 1.15 1.09 0.60 4.37 1.00 1838 1854
## y1[46] -1.15 -1.15 1.16 1.16 -3.05 0.74 1.00 2093 2013
## y1[47] -0.33 -0.33 1.16 1.19 -2.23 1.51 1.00 1935 1888
## y1[48] 2.38 2.39 1.15 1.13 0.51 4.28 1.00 2090 1956
## y1[49] 2.16 2.17 1.12 1.09 0.33 3.98 1.00 2003 2030
## y1[50] 2.33 2.31 1.17 1.12 0.41 4.30 1.00 1715 1714
## m1[1] -0.20 -0.20 0.31 0.31 -0.69 0.31 1.00 1286 1428
## m1[2] -0.26 -0.26 0.26 0.25 -0.72 0.17 1.00 1048 1195
## m1[3] 0.95 0.96 0.29 0.28 0.48 1.44 1.00 1163 1321
## m1[4] 0.19 0.19 0.26 0.25 -0.24 0.63 1.00 1037 1162
## m1[5] -0.06 -0.06 0.32 0.33 -0.59 0.47 1.00 1652 1528
## m1[6] -0.52 -0.52 0.27 0.25 -0.98 -0.08 1.00 1077 1340
## m1[7] 0.52 0.53 0.29 0.30 0.05 0.99 1.00 1566 1557
## m1[8] 2.79 2.79 0.36 0.36 2.18 3.38 1.00 1577 1474
## m1[9] -1.39 -1.40 0.37 0.37 -2.01 -0.79 1.00 1608 1491
## m1[10] -0.26 -0.26 0.33 0.33 -0.81 0.26 1.00 1217 1454
## m1[11] 0.16 0.16 0.26 0.25 -0.27 0.60 1.00 1037 1162
## m1[12] 1.49 1.49 0.26 0.26 1.06 1.91 1.00 1543 1513
## m1[13] 2.43 2.43 0.27 0.28 2.00 2.87 1.00 1780 1521
## m1[14] 1.18 1.18 0.30 0.30 0.69 1.65 1.00 1217 1367
## m1[15] 2.45 2.44 0.35 0.33 1.90 3.03 1.00 1327 1382
## m1[16] 3.42 3.42 0.32 0.33 2.91 3.93 1.00 2143 1484
## m1[17] 2.41 2.42 0.34 0.34 1.85 2.97 1.00 1476 1451
## m1[18] 1.02 1.02 0.30 0.30 0.53 1.50 1.00 1192 1399
## m1[19] -0.10 -0.09 0.33 0.33 -0.64 0.44 1.00 1659 1650
## m1[20] 3.93 3.93 0.36 0.36 3.35 4.51 1.00 2335 1424
## m1[21] -1.94 -1.94 0.34 0.33 -2.49 -1.37 1.00 1487 1550
## m1[22] 0.77 0.77 0.29 0.28 0.30 1.26 1.00 1167 1235
## m1[23] 1.60 1.60 0.32 0.32 1.07 2.12 1.00 1198 1247
## m1[24] 3.17 3.16 0.40 0.39 2.52 3.83 1.00 1463 1526
## m1[25] 1.69 1.69 0.31 0.30 1.20 2.21 1.00 1206 1247
## m1[26] -0.89 -0.90 0.34 0.34 -1.45 -0.33 1.00 1463 1510
## m1[27] -0.46 -0.46 0.27 0.25 -0.91 -0.02 1.00 1068 1251
## m1[28] 2.47 2.47 0.27 0.28 2.03 2.91 1.00 1790 1454
## m1[29] 2.66 2.66 0.28 0.28 2.21 3.12 1.00 1859 1454
## m1[30] 2.29 2.28 0.34 0.32 1.75 2.86 1.00 1296 1370
## m1[31] -1.54 -1.54 0.31 0.30 -2.05 -1.01 1.00 1326 1501
## m1[32] -2.16 -2.15 0.46 0.46 -2.92 -1.43 1.00 1589 1548
## m1[33] 2.15 2.16 0.36 0.35 1.58 2.76 1.00 1310 1346
## m1[34] 0.53 0.53 0.29 0.30 0.05 1.00 1.00 1565 1557
## m1[35] 2.11 2.11 0.33 0.33 1.58 2.63 1.00 1401 1482
## m1[36] -1.08 -1.07 0.40 0.41 -1.74 -0.42 1.00 1812 1608
## m1[37] -0.65 -0.65 0.27 0.25 -1.11 -0.20 1.00 1098 1360
## m1[38] -0.02 -0.02 0.26 0.25 -0.47 0.41 1.00 1037 1162
## m1[39] -2.59 -2.60 0.39 0.39 -3.23 -1.96 1.00 1746 1578
## m1[40] 0.22 0.22 0.29 0.29 -0.26 0.71 1.00 1209 1340
## m1[41] 0.58 0.57 0.27 0.26 0.12 1.02 1.00 1059 1224
## m1[42] -0.64 -0.64 0.37 0.37 -1.25 -0.04 1.00 1756 1660
## m1[43] 2.58 2.58 0.35 0.35 1.99 3.14 1.00 1520 1473
## m1[44] 1.76 1.77 0.26 0.27 1.34 2.19 1.00 1583 1545
## m1[45] 2.47 2.46 0.35 0.34 1.92 3.06 1.00 1331 1382
## m1[46] -1.15 -1.15 0.29 0.28 -1.64 -0.66 1.00 1209 1457
## m1[47] -0.33 -0.34 0.26 0.25 -0.79 0.10 1.00 1054 1226
## m1[48] 2.37 2.37 0.27 0.27 1.94 2.80 1.00 1755 1541
## m1[49] 2.15 2.16 0.26 0.27 1.73 2.58 1.00 1687 1552
## m1[50] 2.38 2.39 0.34 0.34 1.81 2.93 1.00 1467 1388
fn(f2)## Initial log joint probability = -82.744
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 25 -24.0366 0.000100239 0.00122831 1 1 30
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.
## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
##
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -27.29 -26.94 1.91 1.63 -31.04 -24.93 1.01 638 850
## b[1] -0.30 -0.29 0.29 0.28 -0.81 0.18 1.00 600 862
## b[2] 0.88 0.87 0.10 0.10 0.72 1.03 1.00 1814 1588
## b[3] 1.19 1.19 0.49 0.48 0.40 1.98 1.01 473 824
## b[4] 1.49 1.47 0.44 0.43 0.78 2.20 1.01 461 670
## b[5] -1.27 -1.23 0.68 0.68 -2.40 -0.19 1.01 388 457
## s 1.07 1.06 0.12 0.11 0.90 1.28 1.00 1720 1524
## y1[1] 0.17 0.17 1.11 1.10 -1.67 1.99 1.00 2006 1909
## y1[2] -0.49 -0.48 1.12 1.09 -2.33 1.31 1.00 1601 2013
## y1[3] 1.23 1.24 1.12 1.08 -0.56 3.02 1.00 1692 1997
## y1[4] -0.10 -0.10 1.11 1.06 -1.94 1.75 1.00 1870 1973
## y1[5] -0.19 -0.19 1.12 1.05 -2.06 1.70 1.00 2056 2059
## y1[6] -0.77 -0.77 1.11 1.06 -2.60 1.01 1.00 1735 1954
## y1[7] 0.29 0.30 1.13 1.15 -1.56 2.21 1.00 2024 1783
## y1[8] 3.15 3.15 1.13 1.09 1.27 4.93 1.00 1889 1820
## y1[9] -0.97 -0.97 1.18 1.15 -2.91 0.96 1.00 1669 1407
## y1[10] 0.23 0.23 1.16 1.18 -1.68 2.14 1.00 1740 1954
## y1[11] -0.09 -0.11 1.09 1.06 -1.84 1.76 1.00 1735 2001
## y1[12] 1.25 1.27 1.08 1.08 -0.48 3.03 1.00 2032 1876
## y1[13] 2.16 2.14 1.11 1.12 0.35 4.01 1.00 1774 1884
## y1[14] 1.62 1.62 1.14 1.11 -0.15 3.50 1.00 1903 1761
## y1[15] 2.68 2.69 1.16 1.13 0.74 4.62 1.00 1970 1971
## y1[16] 3.10 3.12 1.13 1.13 1.24 4.91 1.00 1862 1739
## y1[17] 2.78 2.78 1.11 1.06 0.93 4.63 1.00 1823 2055
## y1[18] 1.45 1.43 1.15 1.10 -0.38 3.32 1.00 1679 1916
## y1[19] -0.28 -0.27 1.13 1.09 -2.15 1.52 1.00 2051 2127
## y1[20] 3.56 3.56 1.12 1.09 1.69 5.38 1.00 1977 1954
## y1[21] -2.14 -2.16 1.14 1.10 -3.98 -0.26 1.00 1916 1840
## y1[22] 1.08 1.08 1.11 1.09 -0.78 2.90 1.00 1925 1968
## y1[23] 1.26 1.25 1.13 1.13 -0.62 3.10 1.00 1820 1787
## y1[24] 3.39 3.40 1.15 1.12 1.40 5.23 1.00 1962 1876
## y1[25] 1.94 1.96 1.10 1.07 0.12 3.75 1.00 1661 1703
## y1[26] -0.50 -0.50 1.13 1.12 -2.38 1.35 1.00 1711 1786
## y1[27] -0.71 -0.72 1.12 1.11 -2.50 1.16 1.00 1690 1627
## y1[28] 2.19 2.23 1.10 1.06 0.34 3.99 1.00 1977 1895
## y1[29] 2.34 2.37 1.10 1.15 0.54 4.10 1.00 1785 1885
## y1[30] 2.53 2.55 1.10 1.06 0.70 4.32 1.00 1867 1700
## y1[31] -1.73 -1.72 1.11 1.09 -3.59 0.07 1.00 2087 2030
## y1[32] -1.57 -1.60 1.18 1.17 -3.51 0.39 1.00 1847 2009
## y1[33] 1.78 1.82 1.14 1.14 -0.14 3.60 1.00 1836 1745
## y1[34] 0.33 0.35 1.13 1.13 -1.53 2.17 1.00 2164 1930
## y1[35] 2.52 2.55 1.16 1.15 0.64 4.37 1.00 1497 1911
## y1[36] -1.18 -1.20 1.16 1.18 -3.03 0.76 1.00 1988 1894
## y1[37] -0.92 -0.90 1.11 1.14 -2.72 0.91 1.00 1490 1530
## y1[38] -0.30 -0.30 1.15 1.15 -2.12 1.63 1.00 1445 1676
## y1[39] -2.74 -2.73 1.14 1.10 -4.70 -0.85 1.00 1747 1854
## y1[40] 0.59 0.57 1.12 1.09 -1.21 2.40 1.00 1892 2012
## y1[41] 0.27 0.30 1.14 1.08 -1.65 2.12 1.00 1637 1743
## y1[42] -0.76 -0.79 1.11 1.13 -2.57 1.07 1.00 2172 2104
## y1[43] 2.94 2.94 1.11 1.09 1.12 4.78 1.00 2002 1893
## y1[44] 1.51 1.53 1.10 1.08 -0.34 3.35 1.00 2201 2060
## y1[45] 2.72 2.73 1.09 1.09 0.93 4.55 1.00 1819 2050
## y1[46] -1.39 -1.39 1.11 1.09 -3.25 0.46 1.00 2211 1845
## y1[47] -0.62 -0.62 1.11 1.15 -2.42 1.18 1.00 2023 1767
## y1[48] 2.09 2.10 1.13 1.11 0.24 3.94 1.00 1961 2098
## y1[49] 1.87 1.86 1.10 1.07 0.08 3.66 1.00 1938 1919
## y1[50] 2.72 2.72 1.17 1.16 0.81 4.65 1.00 2026 1901
## m1[1] 0.17 0.17 0.34 0.34 -0.40 0.74 1.00 764 1022
## m1[2] -0.52 -0.51 0.29 0.28 -1.02 -0.04 1.00 627 929
## m1[3] 1.26 1.26 0.31 0.30 0.74 1.77 1.00 818 1072
## m1[4] -0.08 -0.07 0.30 0.29 -0.60 0.41 1.00 582 867
## m1[5] -0.22 -0.22 0.33 0.31 -0.76 0.32 1.00 2036 1508
## m1[6] -0.76 -0.76 0.29 0.28 -1.27 -0.29 1.00 666 979
## m1[7] 0.34 0.34 0.30 0.29 -0.16 0.82 1.00 1888 1591
## m1[8] 3.14 3.13 0.38 0.37 2.52 3.78 1.00 1218 1330
## m1[9] -0.97 -0.97 0.41 0.40 -1.63 -0.28 1.00 836 1029
## m1[10] 0.24 0.23 0.39 0.38 -0.38 0.90 1.00 872 1056
## m1[11] -0.12 -0.11 0.30 0.29 -0.63 0.38 1.00 584 868
## m1[12] 1.25 1.26 0.29 0.28 0.77 1.72 1.00 1676 1404
## m1[13] 2.15 2.16 0.31 0.29 1.64 2.64 1.00 1529 1534
## m1[14] 1.61 1.60 0.35 0.34 1.04 2.20 1.00 895 1191
## m1[15] 2.68 2.68 0.35 0.35 2.09 3.24 1.00 1152 1630
## m1[16] 3.08 3.09 0.36 0.34 2.49 3.67 1.00 1457 1490
## m1[17] 2.78 2.77 0.37 0.36 2.20 3.40 1.00 1097 1247
## m1[18] 1.46 1.45 0.36 0.34 0.90 2.05 1.00 881 1226
## m1[19] -0.25 -0.26 0.33 0.31 -0.80 0.29 1.00 2049 1453
## m1[20] 3.57 3.56 0.39 0.38 2.91 4.20 1.00 1454 1394
## m1[21] -2.11 -2.10 0.34 0.33 -2.67 -1.54 1.00 1135 1400
## m1[22] 1.09 1.09 0.31 0.30 0.56 1.61 1.00 794 1068
## m1[23] 1.25 1.27 0.36 0.35 0.65 1.83 1.00 598 868
## m1[24] 3.36 3.36 0.38 0.39 2.72 3.99 1.00 1418 1610
## m1[25] 1.96 1.96 0.32 0.32 1.42 2.48 1.00 941 1319
## m1[26] -0.49 -0.49 0.37 0.37 -1.11 0.13 1.00 797 956
## m1[27] -0.70 -0.69 0.29 0.28 -1.21 -0.22 1.00 656 987
## m1[28] 2.18 2.19 0.31 0.29 1.67 2.67 1.00 1526 1534
## m1[29] 2.36 2.37 0.32 0.30 1.85 2.87 1.00 1511 1484
## m1[30] 2.53 2.53 0.34 0.35 1.95 3.08 1.00 1104 1592
## m1[31] -1.73 -1.72 0.32 0.30 -2.26 -1.21 1.00 959 1187
## m1[32] -1.56 -1.57 0.52 0.51 -2.38 -0.68 1.00 998 1302
## m1[33] 1.78 1.79 0.39 0.40 1.12 2.43 1.00 637 943
## m1[34] 0.34 0.34 0.30 0.29 -0.15 0.83 1.00 1888 1591
## m1[35] 2.50 2.49 0.36 0.35 1.93 3.11 1.00 1025 1296
## m1[36] -1.18 -1.19 0.39 0.38 -1.83 -0.54 1.00 2236 1548
## m1[37] -0.88 -0.88 0.29 0.28 -1.38 -0.41 1.00 690 964
## m1[38] -0.29 -0.28 0.29 0.28 -0.80 0.19 1.00 599 862
## m1[39] -2.73 -2.72 0.37 0.37 -3.35 -2.10 1.00 1431 1464
## m1[40] 0.57 0.57 0.33 0.32 0.02 1.11 1.00 766 1101
## m1[41] 0.28 0.30 0.31 0.30 -0.24 0.78 1.00 568 814
## m1[42] -0.77 -0.78 0.36 0.34 -1.37 -0.17 1.00 2187 1456
## m1[43] 2.94 2.93 0.37 0.36 2.34 3.57 1.00 1165 1276
## m1[44] 1.51 1.53 0.29 0.27 1.03 1.99 1.00 1618 1569
## m1[45] 2.70 2.70 0.35 0.35 2.11 3.26 1.00 1158 1630
## m1[46] -1.36 -1.36 0.30 0.29 -1.88 -0.87 1.00 822 1052
## m1[47] -0.58 -0.58 0.29 0.28 -1.09 -0.11 1.00 637 938
## m1[48] 2.09 2.09 0.31 0.29 1.58 2.57 1.00 1540 1534
## m1[49] 1.88 1.89 0.30 0.28 1.38 2.36 1.00 1568 1612
## m1[50] 2.75 2.74 0.37 0.36 2.17 3.37 1.00 1088 1280
tb=tibble(xa=runif(n,-2,2),xb=runif(n,-2,2),
ca=sample(c('a','b'),n,replace=T),cb=sample(c('a','b'),n,replace=T),
y=rnorm(n,xa+xb-xa*xb+(ca=='b')*2+(cb=='b')*2-(ca=='b')*(cb=='b'),1))
grid.arrange(qplot(data=tb,xa,y,col=ca),
qplot(data=tb,xa,y,col=cb),
qplot(data=tb,xb,y,col=ca),
qplot(data=tb,xb,y,col=cb))
fn=function(f){
X=model.matrix(f,tb)
data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X)
mle=mdl$optimize(data=data)
mle
mcmc=goMCMC(mdl,data)
seeMCMC(mcmc,exc='m',ch=F)
y1=mcmc$draws('y1')
smy=summary(y1)
tb$y1=smy$median
tb$l5=smy$q5
tb$u5=smy$q95
grid.arrange(
qplot(data=tb,y,y1,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))
par(mfrow=c(1,2))
hist(tb$y-tb$y1,xlab='obs.-prd.',main='residual')
density(tb$y-tb$y1) |> plot(xlab='obs.-prd.',main='residual')
tb=arrange(tb,y1)
grid.arrange(
qplot(data=tb,1:n,y,col=ca)+
geom_point(aes(x=1:n,y=y1),tb,col='black')+
geom_line(aes(x=1:n,y=l5),tb,col='gray')+
geom_line(aes(x=1:n,y=u5),tb,col='gray'),
qplot(data=tb,1:n,y,col=cb)+
geom_point(aes(x=1:n,y=y1),tb,col='black')+
geom_line(aes(x=1:n,y=l5),tb,col='gray')+
geom_line(aes(x=1:n,y=u5),tb,col='gray'),
ncol=2
)
}
f0=formula(y~xa+xb+ca+cb)
f1=formula(y~xa+xb+ca*cb)
f2=formula(y~xa*xb+ca*cb)
fn(f0)## Initial log joint probability = -7137.98
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 22 -43.6763 0.00019513 0.00152004 1 1 36
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.
## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.3 seconds.
##
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -46.54 -46.18 1.88 1.68 -50.15 -44.14 1.00 701 926
## b[1] 0.42 0.41 0.38 0.38 -0.19 1.04 1.00 671 1000
## b[2] 0.64 0.64 0.19 0.19 0.32 0.95 1.00 2876 1409
## b[3] 0.70 0.70 0.21 0.19 0.36 1.04 1.00 2291 1320
## b[4] 0.84 0.83 0.45 0.46 0.11 1.58 1.00 772 965
## b[5] 1.94 1.95 0.45 0.43 1.20 2.70 1.01 725 1008
## s 1.58 1.56 0.18 0.17 1.31 1.90 1.00 2042 1585
## y1[1] 2.72 2.71 1.72 1.70 0.00 5.56 1.00 1977 1844
## y1[2] 1.31 1.30 1.63 1.64 -1.34 3.95 1.00 1829 1821
## y1[3] 1.97 1.99 1.71 1.64 -0.90 4.80 1.00 1681 1965
## y1[4] 5.10 5.12 1.72 1.65 2.20 7.84 1.00 1928 1801
## y1[5] 0.15 0.12 1.63 1.60 -2.47 2.85 1.00 1674 1882
## y1[6] 3.91 3.93 1.64 1.61 1.17 6.57 1.00 2045 1969
## y1[7] 3.44 3.44 1.64 1.68 0.73 6.10 1.00 2038 1734
## y1[8] -0.48 -0.46 1.73 1.67 -3.30 2.30 1.00 1988 1945
## y1[9] 4.08 4.05 1.70 1.66 1.29 6.96 1.00 2164 1934
## y1[10] 1.17 1.18 1.66 1.68 -1.53 3.82 1.00 1874 1931
## y1[11] 1.62 1.61 1.64 1.63 -1.14 4.39 1.00 2033 1995
## y1[12] 1.49 1.46 1.66 1.63 -1.33 4.30 1.00 1852 1881
## y1[13] 1.26 1.26 1.64 1.69 -1.43 3.92 1.00 1970 1931
## y1[14] 2.62 2.63 1.69 1.65 -0.18 5.38 1.00 1995 1857
## y1[15] 2.22 2.22 1.70 1.67 -0.52 4.93 1.00 1941 1770
## y1[16] -0.93 -0.97 1.64 1.60 -3.61 1.83 1.00 1858 1940
## y1[17] 0.42 0.45 1.65 1.63 -2.31 3.11 1.00 1953 1896
## y1[18] -1.00 -0.96 1.69 1.63 -3.84 1.66 1.00 2093 1915
## y1[19] 2.30 2.33 1.65 1.65 -0.45 4.96 1.00 1895 1997
## y1[20] 2.40 2.40 1.72 1.65 -0.42 5.24 1.00 2046 1958
## y1[21] 2.83 2.80 1.68 1.62 -0.05 5.57 1.00 1972 1969
## y1[22] 2.73 2.73 1.64 1.65 0.08 5.43 1.00 1957 1919
## y1[23] -0.87 -0.87 1.66 1.68 -3.54 1.84 1.00 1781 1756
## y1[24] 1.69 1.71 1.73 1.70 -1.14 4.62 1.00 1878 1496
## y1[25] 3.87 3.84 1.63 1.60 1.29 6.64 1.00 1822 1972
## y1[26] 2.21 2.17 1.63 1.60 -0.48 4.90 1.00 1814 1742
## y1[27] 1.17 1.16 1.60 1.61 -1.41 3.86 1.00 1980 1879
## y1[28] 3.89 3.88 1.67 1.62 1.10 6.63 1.00 2045 2013
## y1[29] 0.96 0.92 1.71 1.70 -1.82 3.66 1.00 1835 2012
## y1[30] 3.16 3.11 1.61 1.60 0.51 5.88 1.00 2156 2058
## y1[31] 0.61 0.60 1.66 1.60 -2.05 3.34 1.00 2026 1722
## y1[32] 0.90 0.93 1.67 1.65 -1.84 3.62 1.00 1924 1765
## y1[33] -0.46 -0.45 1.65 1.68 -3.16 2.20 1.00 1975 1693
## y1[34] -0.96 -1.00 1.69 1.68 -3.67 1.84 1.00 1825 1835
## y1[35] 3.57 3.59 1.66 1.64 0.88 6.26 1.00 2024 2014
## y1[36] 2.94 2.94 1.68 1.66 0.13 5.73 1.00 2059 1932
## y1[37] 0.59 0.62 1.61 1.57 -2.11 3.22 1.00 1897 1857
## y1[38] 2.02 2.03 1.64 1.53 -0.64 4.74 1.00 1847 1721
## y1[39] 0.66 0.68 1.68 1.67 -2.11 3.36 1.00 1595 1888
## y1[40] 2.03 2.03 1.66 1.64 -0.66 4.56 1.00 1898 1916
## y1[41] 0.14 0.16 1.58 1.61 -2.52 2.62 1.00 1949 1778
## y1[42] 2.12 2.08 1.66 1.68 -0.53 4.80 1.00 1910 1972
## y1[43] 3.32 3.35 1.66 1.57 0.62 6.13 1.00 2055 1714
## y1[44] 3.31 3.35 1.66 1.62 0.59 6.01 1.00 1929 1511
## y1[45] 4.81 4.81 1.67 1.63 2.09 7.51 1.00 2110 1931
## y1[46] 3.90 3.88 1.71 1.64 1.03 6.70 1.00 1860 1749
## y1[47] 2.19 2.24 1.66 1.62 -0.61 4.85 1.00 1660 1902
## y1[48] 1.52 1.54 1.65 1.59 -1.25 4.19 1.00 1973 1786
## y1[49] 1.67 1.68 1.67 1.66 -1.13 4.45 1.00 1802 1894
## y1[50] 2.85 2.88 1.70 1.73 0.05 5.59 1.00 2096 1973
## m1[1] 2.74 2.74 0.61 0.58 1.72 3.74 1.00 1813 1659
## m1[2] 1.33 1.35 0.52 0.51 0.46 2.19 1.00 1493 1162
## m1[3] 1.95 1.95 0.47 0.49 1.20 2.71 1.00 970 1269
## m1[4] 5.09 5.08 0.54 0.56 4.19 5.96 1.00 1645 1280
## m1[5] 0.15 0.15 0.38 0.38 -0.47 0.78 1.00 694 1101
## m1[6] 3.94 3.93 0.52 0.51 3.07 4.81 1.00 1766 1649
## m1[7] 3.44 3.45 0.58 0.57 2.49 4.37 1.00 1780 1519
## m1[8] -0.42 -0.43 0.52 0.51 -1.28 0.47 1.00 1651 1313
## m1[9] 4.05 4.05 0.51 0.51 3.22 4.90 1.00 1703 1600
## m1[10] 1.25 1.26 0.45 0.45 0.47 1.97 1.00 1332 1430
## m1[11] 1.70 1.70 0.47 0.46 0.91 2.43 1.00 1387 1470
## m1[12] 1.50 1.49 0.43 0.43 0.80 2.21 1.00 816 1178
## m1[13] 1.28 1.28 0.45 0.45 0.57 2.03 1.00 1082 1205
## m1[14] 2.64 2.64 0.48 0.46 1.83 3.44 1.00 1495 1622
## m1[15] 2.28 2.28 0.53 0.50 1.40 3.16 1.00 1896 1337
## m1[16] -0.95 -0.95 0.50 0.50 -1.79 -0.14 1.00 1075 1351
## m1[17] 0.41 0.41 0.43 0.42 -0.30 1.14 1.00 1295 1616
## m1[18] -0.99 -0.99 0.49 0.48 -1.79 -0.21 1.00 965 933
## m1[19] 2.27 2.27 0.46 0.44 1.55 3.03 1.00 1314 1297
## m1[20] 2.39 2.40 0.70 0.70 1.21 3.51 1.00 1559 1159
## m1[21] 2.87 2.87 0.51 0.48 2.03 3.71 1.00 1824 1265
## m1[22] 2.76 2.78 0.49 0.49 1.95 3.55 1.00 1479 1285
## m1[23] -0.91 -0.91 0.48 0.49 -1.69 -0.13 1.00 985 1257
## m1[24] 1.67 1.67 0.61 0.57 0.70 2.72 1.00 1433 1169
## m1[25] 3.83 3.83 0.49 0.48 3.06 4.62 1.00 1427 1680
## m1[26] 2.23 2.23 0.50 0.52 1.42 3.04 1.00 1127 1283
## m1[27] 1.17 1.18 0.42 0.41 0.46 1.88 1.00 1268 1124
## m1[28] 3.90 3.91 0.51 0.49 3.04 4.73 1.00 1715 1127
## m1[29] 0.96 0.96 0.57 0.56 0.02 1.91 1.01 924 922
## m1[30] 3.20 3.20 0.49 0.45 2.40 4.00 1.00 1733 1130
## m1[31] 0.53 0.54 0.48 0.44 -0.25 1.32 1.00 1389 1592
## m1[32] 0.88 0.89 0.44 0.42 0.15 1.60 1.00 1288 1595
## m1[33] -0.50 -0.50 0.45 0.45 -1.25 0.26 1.00 969 1414
## m1[34] -0.98 -0.99 0.57 0.56 -1.91 -0.06 1.00 1028 1219
## m1[35] 3.60 3.60 0.65 0.64 2.57 4.67 1.00 1387 1468
## m1[36] 2.92 2.94 0.45 0.42 2.18 3.65 1.00 1574 1131
## m1[37] 0.54 0.53 0.37 0.38 -0.06 1.16 1.00 671 1006
## m1[38] 2.04 2.03 0.42 0.43 1.37 2.73 1.00 1150 1308
## m1[39] 0.69 0.68 0.60 0.58 -0.27 1.70 1.00 1596 1253
## m1[40] 2.01 2.00 0.51 0.49 1.18 2.85 1.00 1783 1219
## m1[41] 0.06 0.07 0.41 0.43 -0.60 0.75 1.00 879 1292
## m1[42] 2.08 2.08 0.43 0.43 1.39 2.77 1.00 1194 1309
## m1[43] 3.39 3.41 0.58 0.58 2.45 4.34 1.00 1800 1478
## m1[44] 3.33 3.34 0.45 0.41 2.60 4.06 1.00 1545 1103
## m1[45] 4.83 4.82 0.61 0.61 3.83 5.82 1.00 1916 1609
## m1[46] 3.89 3.88 0.52 0.51 3.04 4.73 1.00 1255 1424
## m1[47] 2.18 2.18 0.41 0.42 1.51 2.86 1.00 1121 1325
## m1[48] 1.48 1.49 0.46 0.45 0.71 2.24 1.00 1348 952
## m1[49] 1.64 1.64 0.54 0.54 0.76 2.52 1.01 926 886
## m1[50] 2.79 2.80 0.54 0.53 1.90 3.68 1.00 1650 1426
fn(f1)## Initial log joint probability = -3972.6
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 27 -43.0984 0.000110501 0.000793672 1 1 40
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.
## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
##
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -46.49 -46.14 1.94 1.77 -50.18 -44.06 1.00 832 1244
## b[1] 0.26 0.25 0.42 0.42 -0.45 0.97 1.00 856 1184
## b[2] 0.58 0.59 0.19 0.20 0.26 0.89 1.00 2297 1670
## b[3] 0.70 0.70 0.20 0.19 0.39 1.03 1.00 2955 1718
## b[4] 1.22 1.22 0.60 0.58 0.26 2.20 1.01 705 910
## b[5] 2.39 2.38 0.67 0.68 1.30 3.47 1.01 596 790
## b[6] -0.94 -0.96 0.91 0.94 -2.41 0.54 1.00 501 689
## s 1.58 1.56 0.18 0.17 1.33 1.89 1.00 2207 1792
## y1[1] 2.88 2.90 1.71 1.69 0.10 5.72 1.00 2016 2010
## y1[2] 1.68 1.70 1.64 1.62 -0.91 4.34 1.00 1802 1852
## y1[3] 1.71 1.73 1.69 1.64 -1.04 4.34 1.00 1888 1895
## y1[4] 4.77 4.72 1.66 1.70 2.09 7.55 1.00 1536 1727
## y1[5] -0.01 -0.04 1.63 1.57 -2.67 2.67 1.00 2055 1877
## y1[6] 4.10 4.13 1.70 1.65 1.32 6.94 1.00 2040 1931
## y1[7] 3.63 3.65 1.71 1.70 0.77 6.42 1.00 1861 1865
## y1[8] -0.17 -0.14 1.65 1.59 -2.91 2.49 1.00 1795 1701
## y1[9] 4.21 4.16 1.68 1.62 1.46 6.93 1.00 1787 1687
## y1[10] 1.49 1.45 1.64 1.62 -1.23 4.16 1.00 1954 1826
## y1[11] 1.83 1.81 1.66 1.67 -0.87 4.52 1.00 1956 1856
## y1[12] 1.29 1.29 1.68 1.66 -1.52 4.04 1.00 1895 2102
## y1[13] 1.04 1.07 1.73 1.71 -1.81 3.82 1.00 1668 1575
## y1[14] 2.86 2.82 1.67 1.64 0.18 5.68 1.00 1729 1832
## y1[15] 2.11 2.08 1.66 1.66 -0.63 4.77 1.00 2175 1836
## y1[16] -1.09 -1.07 1.69 1.68 -3.89 1.62 1.00 1758 1403
## y1[17] 0.67 0.63 1.64 1.68 -1.89 3.33 1.00 1927 1784
## y1[18] -1.13 -1.12 1.65 1.63 -3.92 1.51 1.00 2003 2009
## y1[19] 2.03 2.06 1.70 1.70 -0.75 4.74 1.00 1987 1917
## y1[20] 2.77 2.73 1.76 1.66 -0.01 5.65 1.00 1970 1915
## y1[21] 2.74 2.68 1.65 1.57 0.16 5.45 1.00 1969 1826
## y1[22] 2.91 2.90 1.71 1.68 0.18 5.64 1.00 1939 1729
## y1[23] -1.04 -1.05 1.64 1.60 -3.74 1.63 1.00 1657 1820
## y1[24] 1.39 1.42 1.73 1.73 -1.50 4.20 1.00 2005 1954
## y1[25] 4.07 4.05 1.67 1.69 1.32 6.76 1.00 1722 1577
## y1[26] 1.94 1.99 1.69 1.64 -0.95 4.78 1.00 1982 1704
## y1[27] 1.42 1.39 1.68 1.63 -1.35 4.33 1.00 2016 1974
## y1[28] 3.69 3.70 1.70 1.62 0.85 6.54 1.00 2003 1706
## y1[29] 0.81 0.82 1.67 1.63 -1.94 3.60 1.00 1836 1527
## y1[30] 2.99 3.00 1.68 1.69 0.26 5.78 1.00 1840 1852
## y1[31] 0.72 0.71 1.63 1.59 -1.95 3.41 1.00 1949 1850
## y1[32] 1.06 1.10 1.65 1.62 -1.73 3.73 1.00 1661 1908
## y1[33] -0.67 -0.68 1.63 1.60 -3.30 2.00 1.00 1981 1547
## y1[34] -1.05 -1.07 1.71 1.69 -3.87 1.83 1.00 1969 2001
## y1[35] 3.17 3.14 1.76 1.73 0.25 6.05 1.00 1923 1782
## y1[36] 2.72 2.71 1.68 1.65 0.00 5.48 1.00 1980 1684
## y1[37] 0.39 0.40 1.66 1.65 -2.39 3.09 1.00 1924 1645
## y1[38] 2.34 2.34 1.69 1.75 -0.38 5.08 1.00 1747 1669
## y1[39] 1.03 1.04 1.69 1.67 -1.69 3.87 1.00 1836 1855
## y1[40] 1.91 1.88 1.68 1.71 -0.85 4.67 1.00 2131 1973
## y1[41] -0.11 -0.11 1.62 1.65 -2.64 2.58 1.00 1933 2049
## y1[42] 2.36 2.31 1.66 1.62 -0.33 5.13 1.00 2027 1976
## y1[43] 3.51 3.53 1.69 1.67 0.65 6.21 1.00 2131 2058
## y1[44] 3.11 3.09 1.67 1.68 0.31 5.82 1.00 1901 1894
## y1[45] 5.05 5.03 1.74 1.69 2.19 7.96 1.00 1921 1776
## y1[46] 3.50 3.50 1.71 1.69 0.67 6.23 1.00 1831 1895
## y1[47] 2.43 2.39 1.70 1.60 -0.38 5.31 1.00 1818 1893
## y1[48] 1.74 1.75 1.65 1.59 -1.04 4.42 1.00 1942 1933
## y1[49] 1.53 1.52 1.63 1.62 -1.12 4.16 1.00 1924 1787
## y1[50] 2.86 2.87 1.69 1.64 0.11 5.66 1.00 1753 2015
## m1[1] 2.91 2.92 0.62 0.58 1.87 3.95 1.00 2064 1603
## m1[2] 1.63 1.63 0.58 0.60 0.70 2.58 1.00 1311 1475
## m1[3] 1.74 1.74 0.52 0.53 0.86 2.60 1.00 860 1350
## m1[4] 4.79 4.79 0.59 0.59 3.80 5.76 1.00 1743 1593
## m1[5] 0.00 0.00 0.42 0.42 -0.70 0.71 1.00 906 1103
## m1[6] 4.13 4.13 0.56 0.53 3.23 5.04 1.00 1670 1740
## m1[7] 3.57 3.58 0.59 0.61 2.58 4.53 1.00 1858 1304
## m1[8] -0.13 -0.13 0.60 0.62 -1.12 0.86 1.00 1243 1447
## m1[9] 4.26 4.25 0.56 0.55 3.38 5.16 1.00 1516 1591
## m1[10] 1.42 1.42 0.49 0.50 0.63 2.22 1.00 1294 1480
## m1[11] 1.85 1.84 0.49 0.49 1.03 2.67 1.00 1383 1398
## m1[12] 1.31 1.31 0.47 0.48 0.53 2.09 1.00 838 1270
## m1[13] 1.04 1.04 0.51 0.51 0.17 1.88 1.00 716 1095
## m1[14] 2.86 2.86 0.53 0.52 1.99 3.74 1.00 1309 1229
## m1[15] 2.12 2.11 0.54 0.54 1.23 3.00 1.00 2431 1565
## m1[16] -1.10 -1.10 0.53 0.51 -1.97 -0.21 1.00 1281 1477
## m1[17] 0.66 0.66 0.50 0.52 -0.15 1.46 1.00 1128 1311
## m1[18] -1.11 -1.11 0.51 0.52 -1.94 -0.27 1.00 1396 1722
## m1[19] 2.02 2.03 0.51 0.51 1.17 2.84 1.00 1746 1426
## m1[20] 2.81 2.77 0.80 0.81 1.55 4.14 1.00 1097 1248
## m1[21] 2.69 2.68 0.52 0.50 1.84 3.55 1.00 2327 1655
## m1[22] 2.95 2.94 0.52 0.52 2.10 3.81 1.00 1516 1397
## m1[23] -1.05 -1.05 0.51 0.50 -1.86 -0.20 1.00 1294 1575
## m1[24] 1.40 1.41 0.65 0.65 0.29 2.44 1.00 1888 1492
## m1[25] 4.09 4.07 0.56 0.56 3.20 5.01 1.00 1221 1464
## m1[26] 2.00 2.01 0.56 0.57 1.03 2.92 1.00 856 1405
## m1[27] 1.44 1.43 0.49 0.51 0.64 2.25 1.00 1133 1533
## m1[28] 3.68 3.68 0.53 0.52 2.82 4.54 1.00 2077 1739
## m1[29] 0.86 0.85 0.56 0.57 -0.07 1.80 1.00 1453 1356
## m1[30] 3.00 3.00 0.50 0.49 2.15 3.84 1.00 2153 1742
## m1[31] 0.73 0.72 0.52 0.54 -0.11 1.57 1.00 1296 1433
## m1[32] 1.07 1.06 0.48 0.49 0.29 1.87 1.00 1208 1514
## m1[33] -0.67 -0.68 0.49 0.47 -1.48 0.18 1.00 1076 1188
## m1[34] -1.03 -1.03 0.58 0.60 -1.97 -0.08 1.00 1731 1731
## m1[35] 3.21 3.19 0.72 0.72 2.04 4.38 1.00 1503 1584
## m1[36] 2.72 2.72 0.48 0.47 1.92 3.52 1.00 1988 1651
## m1[37] 0.38 0.37 0.42 0.42 -0.32 1.08 1.00 833 1175
## m1[38] 2.33 2.33 0.52 0.53 1.47 3.22 1.01 919 1267
## m1[39] 0.99 0.98 0.67 0.64 -0.05 2.13 1.00 1162 1407
## m1[40] 1.83 1.83 0.53 0.53 0.94 2.70 1.00 2273 1568
## m1[41] -0.13 -0.13 0.46 0.45 -0.89 0.67 1.00 891 1034
## m1[42] 2.36 2.36 0.52 0.52 1.52 3.23 1.01 961 1213
## m1[43] 3.51 3.53 0.59 0.61 2.52 4.47 1.00 1841 1422
## m1[44] 3.11 3.12 0.48 0.46 2.32 3.93 1.00 1915 1623
## m1[45] 5.04 5.02 0.65 0.65 4.02 6.07 1.00 1727 1725
## m1[46] 3.54 3.53 0.60 0.60 2.53 4.52 1.00 1393 1529
## m1[47] 2.48 2.47 0.52 0.53 1.64 3.36 1.01 903 1239
## m1[48] 1.75 1.75 0.52 0.53 0.92 2.60 1.00 1224 1395
## m1[49] 1.50 1.49 0.55 0.56 0.58 2.38 1.00 1202 1354
## m1[50] 2.91 2.92 0.55 0.55 1.99 3.82 1.00 1617 1334
fn(f2)## Initial log joint probability = -126.472
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 31 -24.9459 4.85521e-05 0.000907066 0.6786 0.6786 35
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.
## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
##
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -29.35 -29.00 2.14 2.01 -33.42 -26.54 1.00 792 1205
## b[1] 0.51 0.51 0.30 0.28 0.03 1.00 1.01 709 1107
## b[2] 0.92 0.92 0.15 0.15 0.67 1.17 1.00 1834 1212
## b[3] 0.79 0.79 0.14 0.14 0.56 1.03 1.00 3182 1847
## b[4] 1.15 1.15 0.43 0.43 0.46 1.88 1.01 511 843
## b[5] 1.98 1.99 0.46 0.45 1.24 2.73 1.01 452 794
## b[6] -0.89 -0.88 0.14 0.14 -1.12 -0.65 1.00 2410 1681
## b[7] -1.04 -1.01 0.66 0.66 -2.15 -0.02 1.01 394 596
## s 1.11 1.10 0.13 0.12 0.93 1.34 1.00 2082 1563
## y1[1] 5.50 5.48 1.27 1.25 3.51 7.61 1.00 1852 2011
## y1[2] 2.99 3.01 1.23 1.21 0.95 5.03 1.00 2075 1933
## y1[3] 1.33 1.35 1.18 1.14 -0.65 3.25 1.00 1959 1648
## y1[4] 3.50 3.52 1.18 1.16 1.53 5.35 1.00 2097 1857
## y1[5] 0.17 0.18 1.14 1.08 -1.68 2.08 1.00 1678 1735
## y1[6] 3.76 3.76 1.23 1.17 1.66 5.86 1.00 1931 1933
## y1[7] 2.07 2.11 1.21 1.20 0.05 4.03 1.00 2000 2059
## y1[8] -1.90 -1.88 1.22 1.19 -3.94 0.14 1.00 2094 1892
## y1[9] 3.45 3.44 1.17 1.16 1.55 5.40 1.00 2050 2012
## y1[10] 2.27 2.27 1.20 1.23 0.31 4.24 1.00 1767 1917
## y1[11] 2.84 2.87 1.16 1.11 0.94 4.79 1.00 1900 1906
## y1[12] 1.34 1.32 1.18 1.15 -0.59 3.33 1.00 1772 1934
## y1[13] 2.02 2.01 1.19 1.21 0.11 3.95 1.00 1905 1968
## y1[14] 3.74 3.75 1.19 1.18 1.83 5.64 1.00 1950 2014
## y1[15] 1.93 1.92 1.19 1.21 0.00 3.87 1.00 1838 1821
## y1[16] -1.65 -1.66 1.16 1.12 -3.53 0.28 1.00 1740 2148
## y1[17] 0.21 0.18 1.14 1.09 -1.72 2.13 1.00 2043 1940
## y1[18] -2.29 -2.27 1.20 1.15 -4.27 -0.27 1.00 1732 2010
## y1[19] 1.00 1.02 1.19 1.20 -0.91 2.97 1.00 1956 2082
## y1[20] 5.28 5.30 1.31 1.27 3.16 7.47 1.00 1979 1831
## y1[21] 3.28 3.29 1.17 1.14 1.30 5.16 1.00 2040 1918
## y1[22] 2.58 2.62 1.17 1.16 0.66 4.50 1.00 1642 1816
## y1[23] -1.84 -1.84 1.16 1.15 -3.68 0.06 1.00 1878 1947
## y1[24] 0.38 0.39 1.18 1.15 -1.56 2.29 1.00 1816 1973
## y1[25] 3.38 3.40 1.17 1.14 1.36 5.28 1.00 1855 2059
## y1[26] 1.26 1.27 1.19 1.11 -0.71 3.23 1.00 1747 1912
## y1[27] 1.80 1.83 1.19 1.19 -0.24 3.68 1.00 2050 2052
## y1[28] 4.39 4.43 1.21 1.20 2.41 6.33 1.00 1863 1931
## y1[29] 2.24 2.25 1.19 1.17 0.34 4.18 1.00 1813 1844
## y1[30] 3.59 3.60 1.19 1.19 1.67 5.61 1.00 2083 2056
## y1[31] 1.04 1.03 1.17 1.18 -0.88 2.93 1.00 1968 1888
## y1[32] 1.45 1.46 1.19 1.22 -0.52 3.34 1.00 1777 1932
## y1[33] -0.51 -0.50 1.14 1.08 -2.36 1.36 1.00 1815 1972
## y1[34] -1.77 -1.78 1.19 1.15 -3.76 0.16 1.00 1862 1941
## y1[35] 5.45 5.46 1.28 1.28 3.40 7.53 1.00 1951 1892
## y1[36] 2.70 2.69 1.16 1.13 0.83 4.67 1.00 1801 1933
## y1[37] 0.65 0.64 1.17 1.15 -1.27 2.56 1.00 1594 1682
## y1[38] 2.13 2.14 1.22 1.19 0.11 4.11 1.00 1953 2088
## y1[39] -0.33 -0.37 1.20 1.17 -2.32 1.65 1.00 1788 1887
## y1[40] 0.71 0.72 1.21 1.19 -1.26 2.69 1.00 1854 1896
## y1[41] 0.56 0.54 1.19 1.18 -1.41 2.52 1.00 1950 1950
## y1[42] 2.29 2.29 1.18 1.14 0.33 4.25 1.00 1918 2009
## y1[43] 2.25 2.21 1.20 1.15 0.28 4.28 1.00 1832 1892
## y1[44] 3.28 3.27 1.18 1.18 1.34 5.21 1.00 1984 1932
## y1[45] 2.63 2.66 1.30 1.34 0.49 4.69 1.00 1897 1748
## y1[46] 4.09 4.08 1.21 1.16 2.09 6.12 1.00 1670 1831
## y1[47] 2.21 2.21 1.19 1.16 0.28 4.20 1.00 1812 1951
## y1[48] 2.50 2.48 1.20 1.17 0.54 4.48 1.00 2010 1973
## y1[49] 1.96 1.94 1.20 1.17 0.02 4.01 1.00 1899 1916
## y1[50] 2.98 2.99 1.17 1.18 1.03 4.85 1.00 1940 1882
## m1[1] 5.53 5.52 0.60 0.60 4.56 6.51 1.00 2471 1644
## m1[2] 2.97 2.96 0.47 0.45 2.19 3.75 1.00 1371 1617
## m1[3] 1.31 1.30 0.37 0.37 0.69 1.91 1.01 768 1431
## m1[4] 3.52 3.51 0.49 0.48 2.74 4.32 1.00 1684 1726
## m1[5] 0.15 0.15 0.30 0.28 -0.34 0.65 1.01 755 963
## m1[6] 3.76 3.76 0.39 0.37 3.09 4.40 1.00 1536 1573
## m1[7] 2.13 2.13 0.47 0.47 1.33 2.90 1.00 1843 1649
## m1[8] -1.92 -1.92 0.52 0.52 -2.78 -1.06 1.00 1448 1301
## m1[9] 3.46 3.46 0.41 0.40 2.77 4.12 1.00 1578 1654
## m1[10] 2.28 2.28 0.37 0.37 1.68 2.89 1.00 1439 1143
## m1[11] 2.84 2.84 0.38 0.39 2.22 3.47 1.00 1525 1416
## m1[12] 1.32 1.32 0.33 0.33 0.77 1.87 1.01 720 1193
## m1[13] 2.06 2.05 0.39 0.38 1.42 2.71 1.00 724 1347
## m1[14] 3.81 3.81 0.40 0.39 3.11 4.44 1.00 1398 1654
## m1[15] 1.94 1.94 0.38 0.36 1.33 2.55 1.00 2245 1679
## m1[16] -1.65 -1.65 0.38 0.38 -2.28 -1.02 1.00 1276 1522
## m1[17] 0.24 0.25 0.37 0.36 -0.37 0.84 1.00 1111 1115
## m1[18] -2.25 -2.25 0.40 0.40 -2.93 -1.60 1.00 1495 1565
## m1[19] 0.99 1.00 0.38 0.39 0.35 1.58 1.00 1717 1575
## m1[20] 5.29 5.31 0.70 0.71 4.15 6.45 1.00 1281 1264
## m1[21] 3.30 3.29 0.38 0.37 2.66 3.93 1.00 2237 1543
## m1[22] 2.55 2.55 0.37 0.36 1.96 3.16 1.00 1555 1590
## m1[23] -1.85 -1.85 0.38 0.38 -2.48 -1.23 1.00 1315 1580
## m1[24] 0.37 0.38 0.47 0.49 -0.42 1.11 1.00 1819 1575
## m1[25] 3.32 3.32 0.41 0.41 2.64 4.00 1.00 1242 1529
## m1[26] 1.26 1.26 0.41 0.42 0.57 1.93 1.01 807 1516
## m1[27] 1.81 1.81 0.36 0.36 1.22 2.42 1.00 1152 1224
## m1[28] 4.41 4.39 0.40 0.40 3.74 5.07 1.00 2143 1790
## m1[29] 2.23 2.23 0.46 0.45 1.48 2.96 1.00 1297 1170
## m1[30] 3.59 3.58 0.37 0.37 2.97 4.20 1.00 2141 1695
## m1[31] 1.07 1.07 0.38 0.38 0.45 1.70 1.00 1279 955
## m1[32] 1.49 1.49 0.35 0.36 0.90 2.07 1.00 1207 897
## m1[33] -0.49 -0.49 0.35 0.34 -1.05 0.08 1.00 956 1336
## m1[34] -1.79 -1.79 0.43 0.44 -2.51 -1.09 1.00 1579 1471
## m1[35] 5.42 5.42 0.61 0.60 4.44 6.44 1.00 1654 1460
## m1[36] 2.70 2.70 0.33 0.32 2.14 3.24 1.00 1905 1607
## m1[37] 0.67 0.66 0.30 0.28 0.18 1.15 1.01 693 1093
## m1[38] 2.10 2.10 0.37 0.36 1.50 2.69 1.00 844 1261
## m1[39] -0.32 -0.32 0.51 0.51 -1.15 0.52 1.00 1188 1482
## m1[40] 0.74 0.74 0.40 0.40 0.08 1.38 1.00 2184 1654
## m1[41] 0.59 0.59 0.34 0.33 0.03 1.16 1.00 795 1268
## m1[42] 2.29 2.30 0.36 0.36 1.69 2.87 1.00 875 1288
## m1[43] 2.30 2.29 0.46 0.46 1.53 3.04 1.00 1781 1568
## m1[44] 3.32 3.31 0.34 0.33 2.76 3.87 1.00 1867 1598
## m1[45] 2.64 2.65 0.60 0.60 1.66 3.61 1.00 1988 1715
## m1[46] 4.10 4.10 0.43 0.43 3.38 4.80 1.00 1209 1663
## m1[47] 2.26 2.26 0.36 0.35 1.66 2.86 1.00 827 1254
## m1[48] 2.50 2.50 0.39 0.39 1.85 3.14 1.00 1262 1437
## m1[49] 1.95 1.95 0.40 0.39 1.30 2.62 1.00 1040 1273
## m1[50] 3.00 3.00 0.38 0.39 2.38 3.65 1.00 1577 1571
